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EXECUTIVE  SUMMARY 


This  is  the  final  report  of  a  three  year  ARIES  research  program  entitled  "Dynamics  of 
Electronically  Excited  Species  in  Gas  and  Condensed  Phases".  This  research  addressed  the  need  of 
the  Air  Force  High  Energy  Density  Materials  program  to  determine  the  lifetimes  of  potential 
advanced  rocket  propellants. 

The  key  objective  was  to  develop  novel  dynamical  methods  to  be  used  in  elucidating  the 
microscopic  dynamics  controlling  lifetimes  important  in  energy  storage  in  isolated  and  matrix- 
embedded  molecules.  This  objective  was  met  by  (1)  developing  novel  semiclassical  methods  and 
combining  them  with  computer  simulation  technology;  and  (2)  using  these  methods  to 
computationally  treat  the  electronically  inelastic  chemistry  of  light  metastables,  including  helium 
and  hydrogen,  in  the  gas  and  condensed  phase. 

This  report  presents  the  new  technologies  and  illustrative  applications  for  which  the 
relevant  potential  energy  surface  information  existed.  These  applications  successfully  determined 
the  mechanisms  controlling  helium  metastable  lifetimes  in  a  transition  from  the  gas  to  the 
condensed  phase. 

The  methods  validated  by  this  research  can  be  applied  to  problems  where  it  is  feasible  to 
obtain  the  relevant  interaction  potentials  and  couplings.  Different  levels  of  implementation  of  the 
semiclassical  methods  are  considered  in  this  report.  The  levels  vary  from  fully  quantitative  to 
semiquantitative  levels  of  accuracy,  and  depend  on  the  dynamical  problem.  In  practice,  the  most 
relevant  few-body  potential  energy  surface  information  can  be  identified  and  obtained.  The 
methods  developed  here  can  then  be  used  to  determine  the  underlying  complex  dynamical  behavior 
of  energetic  molecules. 
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INTRODUCTION 

High  energy  density  material  (HEDM)  candidates  necessarily  involve  excited  or  energetic 
molecular  states.  Their  metastability  in  the  gas  phase  can  arise  from  spontaneous  radiative  decay, 
predissociation,  or  energy  loss  induced  (via  radiative  and  nonradiative  relaxation,  and  reaction)  in  a 
collision  with  the  ambient  gas  molecules.  Any  HEDM  will  be  eventually  used  in  a  condensed  phase 
environment  for  practical  application.  In  the  matrix  environment,  a  variety  of  time  dependent 
processes  can  occur  that  lead  to  decay  of  the  stored  energy.  (1)  the  excitation  may  relax  by  radiative 
decay,  (2)  it  may  decay  due  to  the  relaxation  induced  by  the  coupling  to  the  lattice  phonons,  (3)  it 
may  relax  into  modes  localized  at  the  site  of  the  HEDM  [such  modes  are  strongly  coupled  to  the 
HEDM  states],  (4)  the  stored  excitation  energy  may  be  resonantly  transferred  to  another  excited 
molecule  [this  is  the  elementary  step  that  leads  to  energy  migration],  (5)  the  excitation  may  be 
nonresonantly  transferred  to  a  chemically  different  matrix  species,  and  (6)  a  chemical  reaction  of 
the  excited  HEDM  may  take  place. 

The  main  goal  of  this  research  program  has  been  to  develop  the  capability  to  identify 
important  quenching  pathways  in  gas  and  condensed  phases  using  theoretical  methods.  In  order  to 
be  able  to  treat  electronically  excited  species,  which  are  a  key  source  of  new  metastable  systems, 
we  have  concentrated  on  developing  dynamical  methodologies  for  electronically  inelastic  collision 
problems.  We  have  addressed  gas  and  condensed  phase  lifetime  problems  with  helium  metastables 
as  prototypes. 

Our  key  accomplishments  in  this  research  program  include: 

(1)  development  and  validation  of  general  and  powerful  semiclassical  methods  for  energetic 

species  and  polyatomic  systems, 

(2)  development  of  reduced  heatbath  models  of  condensed  phase  helium, 

(3)  development  of  models  of  condensed  phase  hydrogen  matrix,  and 

(4)  development  of  simulation  procedures  for  solution-phase  reaction  and  cluster  formation 

of  HEDM  with  solvents. 

As  an  application,  we  present  elucidation  of  the  important  dynamical  mechanisms 
contributing  to  the  observed  helium  metastable  atom  quenching  propensities  in  gas  and  condensed 
phases  using  accurate  ab  initio  input  for  the  appropriate  electronic  potentials  and  couplings.  A 
variety  of  processes  have  been  addressed  in  validations  or  illustrative  examples:  (1)  gas  and 
condensed  phase  collision-induced  radiative  quenching,  (2)  gas  and  condensed  phase  collisional 
(nonradiative)  quenching,  (3)  gas  phase  vibrational  energy  transfer,  (4)  gas  phase  electronic  energy 
transfer,  (5)  condensed  phase  reaction  with  solvent  medium,  (6)  pressure  dependent  solvation 
effects  on  different  electronic  states. 
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An  illustrative  example  of  our  results  concerns  the  dynamical  mechanisms  for  the  observed 
lifetime  trends  in  helium  metastables.  Table  1  shows  the  known  lifetimes  of  electronically  excited 
He*(23S)  atoms  and  He2*(a3£u+)  molecules  in  vacuum  and  condensed  phase  environments.  The 

key  trend  to  notice  is  the  large  change  of  lifetime  (from  8000  sec  to  15  (isec)  that  occurs  upon 
transferring  an  excited  atom  into  liquid,  compared  to  the  relatively  small  change  (from  18  sec  to  10 
sec)  found  for  the  diatom.  The  large  change  in  the  vacuum  lifetime  of  the  atomic  versus  diatomic 
species  is  also  noteworthy  since  this  change  is  mostly  accounted  for  by  the  spin- forbidden  radiative 
quenching  process  constantly  operative  in  the  molecule.  Unravelling  the  relevant  mechanisms  to 

Table  1.  Helium  Metastable  Lifetimes 


Type 

Species 

Lifetime 

Basis 

Vacuum 

He*(23S) 

8000  s 

PRL  30,  775  (1973) 

Vacuum 

He2*(a3Xu+) 

18  s 

JCP  90,  2504  (1989) 

Intrinsic  in  4He  Liquid 

He*(23S) 

15  (is  (T-indep) 

PRL  28,  792  (1972) 

Intrinsic  in  4He  Liquid 

He2*(a3V) 

>10  s 

JLTP  36,  47  (1979) 

In  4He  Liquid 

He2*(a3Eu+) 

30  ms 

PRL  28,  792  (1972) 

at  p=1012cnr3 

Timescales 

Key  Process 

Timescale 

He*(23S)  to  He(l’S)  excitation  transfer 

10-8  s 

Cavity  Formation  in  Helium  Liquid 

U5 

i 

o 
**— < 

explain  the  lifetime  changes  of  the  atomic  species  involves  exploring  systematically  the 
electronically  nonadiabatic  dynamical  processes  that  result  in  the  quenching  of  He(23S)  to  He(’S) 
at  the  gaseous  helium  densities  and  at  densities  characteristic  of  condensed  phase  helium 
environments. 
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The  8000  sec  isolated  gas-phase  lifetime  of  He(23S)  is  controlled  by  the  dynamics  of 
atomic  spontaneous  emission  whose  timescale  is  determined  by  relativistically  induced  magnetic 
multipole  couplings  in  the  atom.  The  shift  to  lower  values  of  lifetime  involves  quenching  dynamics 
in  collisions  with  other  helium  atoms.  This  introduces  nonadiabatic  dynamics  involving  scattering 
on  multiple  electronic  surfaces,  notably  the  a(3Eu+),  A(1ZU+),  b(3rig),  B^IIg),  c(3Zg+),  and 

C(1£g+)  states  of  He2.  Dynamical  studies  of  the  important  pathways  leading  to  such  an  overall 

spin-forbidden  quenching  can  now  be  studied  employing  Yarkony's  new  ab  initio  results 
(obtained  in  a  parallel  HEDM  effort)  for  the  relevant  potential  energy  surfaces  and  nonadiabatic 
couplings  of  this  metastable  system.  The  possible  dynamical  pathways  can  be  examined  by 
considering  emission  of  each  of  the  photons  at  many  possible  energies  (frequencies),  thereby 
providing  a  thorough  computational  exploration  of  the  radiative  quenching  spectra.  However,  an 
analysis  of  the  important  pathways  may  simply  be  based  on  selective  dynamical  results  for  key 
scattering  propensities;  in  this  report  we  illustrate  the  usefulness  of  such  selected  studies. 

Radiative  quenching  of  the  gas  phase  He2  (a3lu+)  state  to  the  ground  electronic  state  is  a 

dissociative  process  and  can  be  viewed  as  a  half-collisional  process.  The  gas  phase  lifetime  of  18 
sec  has  been  explained  in  terms  of  radiative  quenching  by  Chabalowski  et  al.  [JCP  90,  2504 
(1989)].  Collisional  quenching  of  the  diatomic  HEDM  by  a  third  He  atom  can  also  be  studied  by 
our  validated  methods.  However,  the  potential  surfaces  and  electronic  couplings  are  not  yet 
available.  Nevertheless,  a  reasonable  understanding  of  what  happens  to  helium  metastables  in 
condensed  phases  emerges  from  our  approximate  pairwise  additive  potential  surfaces  in  the 
computer  simulation  studies  described  below.  This  is  because  helium  metastables  form  bubbles  in 
the  condensed  medium  which  reduce  the  proximity  between  excited  species  and  solvent.  As  seen 
from  Table  1,  the  timescale  for  He*(23S)  to  He(l1S)  excitation  transfer  is  orders  of  magnitude 
larger  than  the  timescale  for  cavity  formation  in  helium  liquid.  This  is  the  reason  why  cavity 
formation  occurs  rather  than  excitation  transfer,  which  would  have  led  to  energy  migration.  The 
formation  of  bubbles  essentially  traps  the  excitation  energy  (thus  making  the  study  of  the  bubble 
structure  and  dynamics  relevant  to  HEDM's). 

The  following  section  of  technical  discussions  is  divided  into  four  main  parts.  Each  part 
contains  the  corresponding  bibliography  at  its  end.  In  Gas  Phase  Methods  and  Calculations, 
formal  developments  involving  semiclassical  computational  methods  are  described,  He-He* 
scattering  calculations  that  include  electronic  inelasticity  are  presented,  and  two  studies  that  involve 
vibrationally  and  vibronically  inelastic  collisions  are  described.  In  the  Condensed  Phase  Modelling 
and  Computer  Experiments  section  computer  simulations  of  the  condensed  phase  structure 
surrounding  metastable  helium  atoms  in  high  pressure  helium  liquid  are  described  and  the 
construction  of  parameters  for  stochastic  heatbath  models  from  such  simulations  illustrated.  This  is 
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followed  by  the  Condensed  Phase  Dynamics  section,  where  the  results  presented  are  based  on 
stochastic  dynamics  calculations  of  helium  metastable  atomic  dynamics  in  a  high-pressure  helium 
matrix  which  forms  microscopic  bubbles  around  the  metastable  excited  atom.  The  fourth  part. 
Extended  Studies,  contains  calculations  of  solvent  shifts  of  the  absorption  line  for  the 
He(23S)->He*(23P)  transition  obtained  from  simulations;  this  part  also  contains  preliminary 
results  from  liquid  hydrogen  simulations. 

TECHNICAL  DISCUSSION 

Gas  Phase  Methods  and  Calculations 

Introduction 

There  is  a  key  need  in  HEDM  research  for  theoretical  methods  of  treating  electronically 
excited  energetic  species  in  many-atom  environments.  Our  efforts  have  developed  novel 
semiclassical  methodologies  that  show  promise  for  use  in  these  applications  (unlike  exact  quantal 
methods,  which  are  prohibitively  expensive  and  usually  intractable  for  problems  of  interest  here). 
We  have  achieved  a  significant  breakthrough  in  the  semiclassical  methodology  during  this  period. 
The  following  subsection  describes  a  new  extended  semiclassical  formalism  that  has  evolved  from 
the  present  computational  research.  Applications  are  then  presented  of  electronically  inelastic 
helium  scattering  problems  in  early  model  studies  as  well  as  studies  employing  the  accurate 
potentials  and  electronic  couplings  that  became  recently  available  for  nonadiabatic  He*-He 
collisions.  This  is  followed  by  a  subsection  on  numerical  benchmark  studies  validating  an 
extremely  efficient  implementation  of  the  present  semiclassical  technology  for  vibrational  excitation 
problems  (in  O+HF  and  Na+H2  collisions).  The  last  subsection  presents  our  quantum  scattering 
work  employed  for  benchmarking  the  Na+H2  system. 

Semiclassical  Methodology  for  Computing  Multichannel  Eikonal 

Wavefunctions  in  Molecular  Collisions:  A  Reformulation  and  Extension 

1.  Introduction 

The  standard  quantum-mechanical  description  of  collisional  processes  is  based  on  the 
adiabatic  separation  of  nuclear  and  electronic  motions,  and  assumes  that  the  electrons 
instantaneously  adjust  to  nuclear  motions  because  of  a  large  difference  in  their  masses.  Encounters 
taking  place  on  a  single  adiabatic  surface  can  be  described  in  the  zeroth -order  approximation  simply 
by  the  classical  equations  of  motion.  However,  transitions  from  one  electronic  surface  to  another 
are  essentially  quantum-mechanical  processes  which  have  no  analogs  in  classical  mechanics.  As 
numerical  solution  of  the  Schroedinger  equation  is  near  the  limit  of  capability  of  modem  computers 
(even  for  four-atom  nonreactive  collisions  on  a  single  potential  surface1)  there  is  presently  no  hope 
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to  obtain  necessary  information  about  electronically  inelastic  scattering  processes  in  polyatomic 
systems  by  performing  exact  quantum-mechanical  calculations. 

The  gap  between  the  quantum  and  classical  treatments  can  however  be  filled  by 
semiclassical  theory.  Construction  of  semiclassical  solutions  is  simpler  compared  with  solving  the 
Schroedinger  equation  because  they  are  often  governed  by  a  manageable  number  of  ordinary 
differential  equations  of  the  first  order  [instead  of  the  (usually  large)  complex  set  of  partial 
differential  equations  exploited  in  quantum  mechanics].  The  semiclassical  theory  is  often 
sufficiently  accurate  to  reproduce  specific  quantum-mechanical  features  which  are  completely 
absent  in  a  purely  classical  picture.  In  particular  the  semiclassical  theory  is  powerful  enough  to 
describe  the  nonadiabatic  transitions  which  are  the  main  concern  of  the  present  research. 

However,  there  is  an  obstacle  that  has  been  damping  progress  in  this  area  for  many  years. 
The  use  of  semiclassical  equations  implies  that  higher-order  corrections  are  negligible  compared 
with  the  kinetic  energy  of  the  system,  so  that  they  break  down  each  time  the  kinetic  energy 
becomes  small,  which  occurs  at  classical  turning  points.  As  a  result,  the  semiclassical 
wavefunctions  become  singular  in  such  regions  The  one-dimensional  problem  has  a  well-known 
solution  within  the  WKB  method2  that  gives  the  recipe  for  continuing  the  semiclassical 
wavefunction  after  its  reflection  at  the  turning  point.  This  recipe,  however,  essentially  exploits  the 
conservation  of  flux  in  connecting  the  incoming  and  outgoing  waves  and  this  is  insufficient  in  the 
multi-dimensional  case  because  of  a  possible  exchange  of  the  probability  density  between  different 
degrees  of  freedom.  The  current  state-of-the-art  theory  in  this  area  is  mainly  developed  at  the  level 
of  formal  theorems3  which  prove  existence  of  the  necessary  solutions  of  the  time-independent 
Schroedinger  equation,  but  do  not  provide  any  practical  numerical  algorithms  to  find  them.  We  can 
cite  only  two  works4’5  where  semiclassical  scattering  wavefunctions  have  been  explicitly 
constructed  in  more  than  one  dimension.  These  works  exploit  completely  different  approaches:  a 
rather  sophisticated  numerical  implementation  of  the  formal  theory  as  it  is  presented  in  Ref.3,  and 
development  of  some  approximate  schemes  still  accurate  enough  to  give  the  necessary  quantitative 
answers.  The  latter  more  pragmatic  view  of  the  theory  presently  seems  to  be  the  only  realistic  way 
to  approach  the  problems  of  practical  interest  Our  present  research  adopts  the  latter  view. 

Until  now,  much  more  attention  has  been  to  given  to  an  alternative  direction  exploiting 
angle-action  variables,  following  the  pioneering  works  of  Marcus6  and  Miller7.  However,  the 
reported  results  in  this  direction  are  very  limited.  Even  disregarding  the  fact  that  the  latter  approach 
is  related  with  the  Schroedinger  equation  only  through  the  correspondence  principle  (and  hence  it 
remains  unclear  to  what  extent  its  predictions  are  equivalent  to  those  of  the  asymptotic  semiclassical 
theory3’8)  there  is  a  purely  practical  reason  for  development  of  the  semiclassical  theory  in  terms  of 
geometrical  variables  --  representation  of  realistic  potential  surfaces  in  angle-action  variables  is  a 
complicated  computational  problem  that  makes  the  whole  scheme  infeasible  for  applications. 
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An  additional  problem  arises  when  the  semiclassical  technique  is  applied  to  nonadiabadc 
processes.  The  problem  arises  from  quantal  interference  between  wavefunctions  propagated  on 
different  potential  surfaces. 

The  present  subsection  reports  some  new  developments  in  overcoming  both  problems  of 
the  nonadiabadc  semiclassical  theory  which  make  its  implementation  possible  for  realistic  systems. 
Assumpting  adiabatic  separation  of  nuclear  motions  in  the  semiclassical  picture  allows  us  to 
connect  incoming  and  outgoing  semiclassical  solutions,  correctly  reproducing  interference  effects 
between  different  trajectories.  By  including  this  adiabatic  hypothesis  into  the  Self-Consistent 
Eikonal  Method  (SCEM)9  we  have  obtained  a  practical  tool  for  treating  electronically  inelastic 
processes  with  a  few  degrees  of  freedom.  It  is  shown  that  SCEM  corresponds  to  a  physical 
picture  in  which  one  approximates  the  flows  of  probability  density  on  different  adiabatic  potential 
surfaces  by  the  average  flow  obtained  by  weighting  the  individulal  flows  with  the  probabilities  for 
the  system  to  be  in  the  appropriate  electronic  state.  We  also  discuss  some  additional  terms  that  were 
previously  ignored  which  arise  in  the  forces  governing  nuclear  motions  for  a  system  with  two  or 
more  degrees  of  freedom  and  suggest  a  numerical  algorithm  for  their  evaluation. 

For  a  more  accurate  treatment  of  electronic  inelasticity  in  molecular  collisions  we  developed 
a  new  extended  approach  called  the  Adiabatic  Velocity  Field  Method  (AVFM).  This  approach  goes 
beyond  the  Ehienfest  effective  potential  approximation  employed  in  SCEM  without  compromising 
on  the  advantages  of  the  latter.  This  extended  theory  also  exploits  the  eikonal  ansatz  for  nuclear 
wavefunctions  in  each  electronic  state,  but  each  wavefunction  is  propagated  now  along  its  own 
trajectory  run  on  the  appropriate  adiabatic  potential  surface.  One  of  the  advantages  of  propagating 
classical  trajectories  on  the  particular  adiabatic  surfaces  as  prescribed  by  AVFM  is  that 
approximations  (in  addition  to  the  short  wavelength  approximation)  are  made  only  in  the 
nonadiabatic  region.  The  equations  outside  the  nonadiabatic  region  are  the  exact  limit  of  the 
Schroedinger  equation  in  the  adiabatic  representation  as  h  tends  to  zero;  hence  their  solutions 
asymptotically  reproduce  specific  features  of  the  exact  wavefunctions  caused  by  local  topology  of 
the  adiabatic  potential  surfaces.  In  particular,  the  new  technique  makes  it  possible  to  treat  tunneling 
along  the  streamline  coordinate  R  in  electronically  inelastic  processes  provided  it  takes  place 
beyond  the  region  of  strong  nonadiabatic  couplings. 

Semiclassical  multichannel-approximation  in  the  mixed  adiabatic/diabatic  representation 

Let  us  consider  the  multichannel  problem  described  by  a  set  of  linear  partial  differential 
equations  of  the  second  order 


14 


(1) 


Y  h2  ]  2- 

-  -y  dQ  -  EJ1  +  V(Q)  -  n  F(Q)  .  V  +  H(Q) 


v  =  o , 


where  V(Q)  is  the  column  formed  by  the  functions  X|/^(Q),  ^(Q)*  —*Vn(Q)  sought  for, 


42=  Ig 


-1/2  JL 


3Q 


g1/2  g 


w*‘  a 


9Q 


(2) 


and  ^  arc  respectively  the  covariant  Laplacian  (see  next  section  for  comments)  and  the  gradient  in 
the  space  of  curvilinear  coordinates  Q.  V(Q)  is  a  n  x  n  diagonal  matrix  having  some  potentials 

Vv(Q)  as  its  diagonal  elements  and  H(Q)  is  the  matrix  of  potential  couplings  with  zeros  at  its  main 

diagonal.  The  components  of  the  vector  F  are  n  x  n  real  antisymmetric  matrices  describing 

nonadiabatic  couplings.  The  set  of  curvilinear  coordinates  Q  is  composed  of  a  streamline 
ft  112  2 

coordinate  Q  =R  and  others  [Q  =w,  Q  sw  ,  ]  describing  some  quasiperiodic  degrees  of 

freedom. 

In  the  adiabatic  representation  Vv(Q)  are  the  adiabatic  potentials  and  the  matrix  H(Q)  =  0 

so  that 
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1  +  V(Q)  -h  F(Q)  •  V 
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The  main  advantage  of  the  adiabatic  representation  is  that  the  last  term  in  brackets  in  Eq.  (3) 
vanishes  in  the  limit  as  h  -»  0.  As  a  result  we  come  to  the  Hamilton-Jacobi  equations  governed  by 
the  adiabatic  potentials.  We  can  thus  treat  the  coupling  term  as  a  perturbation.  We  show  below  that 
one  needs  an  assumption  of  such  a  kind  when  continuing  the  semiclassical  wavefunction  through 
turning  points  of  quasiperiodic  motions.  Another  alternative  is  a  system  with  only  potential 
couplings  which  are  relatively  small  compared  with  the  potentials  Vv(Q).  In  the  one-dimensional 

case  it  is  sufficient  to  require  that  the  couplings  are  negligible  near  the  turning  point  of  the 
streamline  motion.  We  come  to  the  problem  of  this  type,  when  treating  quasiperiodic  degrees  of 
freedom  quantum-mechanically  by  representing  the  Schroedinger  equation 
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in  the  matrix  form 
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assuming  that  the  streamline  coordinate  R  is  orthogonal  to  all  others.  As  a  result  Eq.(4)  takes  the 
form 
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Applying,  by  analogy,  a  similar  expansion  to  an  electronically  inelastic  process  described  by  the 
Schroedinger  equation,  Eq.  (2),  we  come  to  the  set  of  ordinary  differential  equations  in  the  mixed 
representation 
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There  are  two  questions  which  we  need  to  address  here.  First  of  all,  we  construct 
asymptotic  solutions  of  Eq.  (1)  as  h  — >  0-  This  is  a  relatively  easy  problem  which  can  be 
formulated  in  an  arbitrary  set  of  variables.  It  is  much  more  difficult  to  find  the  connection  formulas 
between  different  asymptotic  solutions.  (As  stressed  by  Schiff10  the  asymptotic  solutions  are  of 
little  use  to  us  unless  we  know  how  to  connect  them  together.)  This  part  of  the  problem  can  be 
solved  only  with  some  additional  assumptions  concerning  the  motions  under  discussion.  Here  we 
restrict  our  discussion  to  systems  with  a  single  unbounded  degree  of  freedom  (a  streamline  motion) 
coupled  with  several  quasi-periodic  degrees  of  freedom  such  that  each  motion  has  a  different  time 
scale.  There  may  be  also  some  cyclic  motions  which  do  not  create  any  difficulties  on  their  own.  All 
the  bound  motions  are  uncoupled  at  large  values  of  the  streamline  coordinate  R  and  are  assumed  to 
be  adiabatically  separable  in  the  interaction  region  so  that  one  can  apply  the  usual  WKB 
quantization  rule.2 
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Let  us  illustrate  these  assumptions,  using  scattering  of  a  diatomic  molecule  by  an  infinitely 
massive  uncorrugated  surface  as  an  example.  The  cyclic  motions  are  represented  by  the  precession 
of  the  diatomic  around  the  normal  drawn  from  the  center  of  mass  of  the  diatomic  to  the  surface. 
Relative  oscillations  of  two  atoms  and  the  nutation  of  the  diatom  give  us  two  quasi-periodic 
degrees  of  freedom.  Since  stretching  vibration,  described  by  the  variable  r,  usually  has  a  much 
larger  frequency  than  the  nutation  we  treat  the  former  as  a  quasi-periodic  motion  with  a  frequency 
parametrically  dependent  on  the  nutational  angle  0  and  the  streamline  coordinate  R.  On  the  other 
hand,  the  nutation  is  expected  to  be  faster  than  streamline  motion  so  that  we  can  describe  it  as 
quasi-periodic  motion  with  the  frequency  depending  parametrically  on  R.  We  thus  neglect  the 
effect  of  small-amplitude  oscillations  of  the  interatomic  distance  r  on  the  nutation  of  the  diatomic 
and  treat  the  diatomic  as  a  rigid  linear  rotator  in  this  particular  context.  We  shall  see  below  that  one 
needs  these  assumptions  only  to  formulate  the  connection  formulas  for  the  asymptotic  solutions 
and  we  do  not  make  any  approximations  in  the  classical  Hamiltonian  itself  to  be  consistent  with 
these  assumptions. 

Note  that  we  included  precession  in  our  analysis  only  to  be  able  to  run  trajectories  in  the 
Cartesian  coordinates.  One  can  directly  start  from  the  Schroedinger  equation  expressed  in  terms  of 
the  curvilinear  coordinates  R,  r,  0.  The  appropriate  kinetic-energy  operator 

X-  h2  32  ft2  3  ,23  h2  1  3  ■  »3 
2H.3R2  "  2nrr2 


is  very  similar  to  that  for  the  J=0  atom-diatom  nonreactive  scattering  problem: 
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where  r  is  the  interatomic  distance  in  the  diatomic  molecule,  y  is  the  angle  between  the  diatomic 
molecule  and  the  vector  R  drawn  from  the  incident  atom  to  the  center  of  mass  of  the  diatomic 
molecule,  and  R=l  R I.  (Note  that  the  volume  elements  dR  dr  sin0  d0  and  dR  dr  siny  dy  are  defined 

in  both  cases  in  exactly  the  same  way.)  Therefore  vibrational  excitation  of  a  diatomic  molecule  in  a 
J=0  collision  with  an  atom  and  scattering  of  diatomic  molecules  by  incorrugated  surfaces  can  be 
formally  treated  in  terms  of  the  nearly  identical  formalisms. 
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To  extend  the  theory  of  semiclassical  transition  probabilities  to  the  multichannel  equation, 
Eq.  (1),  we  start  from  Dirac-Marcus'  representation6’11  of  the  wave  functions  \|/j(Q),  ^(Q),  — 

Vn(Q) as 

¥V(Q)  =  Av(Q)exp[iwv(Q)/h]  ,  v  =  l,2,...,n  (10) 

where  A  (Q)  and  W  y  (Q)  are  real  functions  which  vary  slowly  with  their  arguments. 

Substituting  (2.10)  into  (2.1)  with 
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and  neglecting  the  terms  of  the  order  of  h  ,  we  come  to  the  following  set  of  differential  equations 
of  the  first  order. 
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referred  to  below  as  the  Semiclassical  Multichannel  Equations  (SMEs).  Note  that  div  here 

— > 

implies  use  of  the  covariant  derivatives  of  a  vector  I  so  that 
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One  can  easily  verify  that 
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and  hence  the  function  W  (Q)  is  Hamilton's  characteristic  function12  for  the  effective  potential, 
given  by  Eq.  (15). 

The  SMEs  have  the  trivial  solution  for  the  single-channel  problem: 


+  V1(Q)  =  E  , 


(18a) 


div 


with 


P  j(Q)  =  Aj(Q)  . 


(18b) 


(19) 


However,  what  is  taken  for  free  in  that  case  becomes  a  rather  complicated  problem  if  one 
has  to  include  interference  between  the  channels.  Even  the  initial  incoming  values  of  effective 
potentials  Eq.  (16)  cannot  be  unambiguously  determined  because  the  boundary  condition  for 
Hamilton's  principal  function  W  v  (Q)  is  known  only  for  the  initially  populated  channel  labeled 

below  by  index  1.  As  a  result,  quantum-mechanical  correction  Eq.  (17)  to  the  potential  governing 
classical  trajectories  turns  out  to  be  an  ill-defined  function  in  the  infinite-separation  limit 

One  of  the  ways  to  bypass  this  difficulty  is  to  neglect  quantum-mechanical  correction  Eq. 
(17),  compared  with  the  adiabatic  potential  Vv(Q),  bearing  in  mind  that  this  correction  is 
proportional  to  h  and  hence  it  disappears  in  the  limit  h  — >  0 .  We  refer  to  this  approximation  as  the 
Adiabatic  Velocity  Field  Method  (AVFM).  In  the  one-dimensional  case  the  neglect  of  correction 
Eq.  (17)  leads  us  to  the  set  of  ordinary  differential  equations 
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with 


P  (R)=  /2[E— Vv(R)] 


vL 


that  can  be  solved  by  means  of  the  finite-difference  method  provided  that  the  coupling  coefficients 
are  negligibly  small  near  the  turning  points  and  hence  one  can  make  use  of  the  standard  WKB 
connection  formulae  in  each  channel.  The  same  assumption  is  used  by  us  for  the  streamline  motion 
in  the  multi-dimensional  case;  however  the  necessity  to  solve  the  set  of  partial  differential  equations 
makes  the  problem  much  more  complicated. 


Let  us  now  express  each  of  Eqs.  (13)  in  terms  of  its  own  set  of  the  curvilinear  coordinates 
V)  with  q^  (s>0)  used  for  the  initial  values  w^  of  the  quasi-periodic  coordinate  ws  on 

trajectories  governed  by  the  potential  V  (Q)  and  parametrized  by  the  time  t(  v )  =  q® different 

for  each  channel.  Making  use  of  the  well-known  expression  for  the  covariant  derivative  (see 
Eq.(6.87)  of  Kyrala13  with  g1/2  for  J  here): 
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where  g(v)  are  the  coefficients  of  momentum  coupling  in  the  new  set  of  coordinates,  1  v  is  the 
Jacobian  of  the  transformation  from  q^  v )  to  Q  ('positive  bv  definition!  and 
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Taking  into  account  that 
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we  can  represent  Eq.  (13)  as 
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(28) 


(29) 


with  index  1  used  for  the  initially  populated  channel.  Substituting  Eqs.  (16),  (17)  and  (29)  in  (28) 
we  find  the  following  equation  for  the  phase  Wy: 


§=2[E-vf(Q(v)[.(vyWo,)] 

and  hence 

*(v) 

wv»(vrwo>-2  f  dt'{E-vvff(Q(vft(vr»oj} 


(30) 


(31) 


21 


or,  in  a  more  familiar  form: 


t(v)  (v)  U 

Wv(t(v)’W0)=  dQ(v) 


(32) 


Note  that  SMEs  (13)  are  independent  of  the  particular  choice  of  the  curvilinear  coordinates 
Q  as  we  assumed  that  the  semiclassical  approximation  can  be  directly  applied  to  covariant  Laplacian 
Eq.  (2).  Use  of  the  Podolsky  transformation14  to  get  rid  of  the  weight  function  in  the  volume 
element  (represented  in  terms  of  the  curvilinear  coordinates  in  question)  leads  one  to  a  quantum- 
mechanical  correction  of  the  order  of  h  which  does  depend  on  the  particular  choice  of  variables. 
The  correction  is  negligible  compared  with  the  adiabatic  potential  unless  it  is  singular  in  the 
classically  allowed  region.  For  problems  of  interest  here  such  complications  come  from  nutational 
motion.  Namely,  if  the  Podolsky  transformation  is  applied  to  the  kinetic  energy  operator  Eq.  (8), 
this  correction  has  the  form 


5U  =  - 


1 

sin20 


(33) 


As  it  becomes  negatively  infinite  at  0  =  0,ji  the  effective  potential  in  absence  of  precession  has 
infinitely  deep  holes  which  may  trap  trajectories.  When  discussing  free  nutational  motion  Landau 
and  Lifshitz  merely  neglect  this  term  by  narrowing  the  region  of  validity  of  the  semiclassical 
approximation.  There  is  however  a  much  more  ponderable  argument  against  including  this  term  in 
the  effective  potential.  As  primarily  stressed  by  Langer15,  to  apply  the  WKB  method  one  should 
first  express  the  Schroedinger  equation  in  terms  of  a  variable  which  changes  from  -  oo  to  +  oo .  We 
show  in  Appendix  A  that  no  singular  term  appears  if  the  interval  (0,  tt )  is  mapped  onto  the  x-axis 
by  the  transformation  x  =  In  tan  (0/2) ,  and  that  the  Bohr-Sommerfeld  quantization  rule  gives  the 
exact  result  for  rotational  frequencies  when  applied  to  free  nutational  motion  described  in  terms  of 
the  new  variable.  The  direct  consequence  of  our  analysis  is  that  the  kinetic  energy  operator  for  a 
diatom  scattered  by  a  incorrugated  surface  should  be  taken  in  the  form: 

j  =  _  ft  3 _ ft  3  r?.B  h  3  ,  h  m^ 

2^r  9R2  2prr2  to  to  2nrr2  ae2  2nrr2sin20 
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whereas  for  the  J=0  atom-diatom  nonreactive  scattering  problem: 


T  = - ^r4rR 


2  9  h2  a  Tld  h 2(  1 

2^  R2  9R-  9R-24rr2ar  3r  "  2 


1  92 


4dr  4rr 


dy 


2  ‘ 


(35) 


Rotational  energies  for  initial  and  final  states  are  given  by  the  semiclassical  expression  Eq.  (94). 

To  evaluate  the  Jacobian  Jy  we  make  use  of  the  numerical  algorithm  developed  by  Stodden 

and  Micha16,  extending  it  to  the  equations  of  motion  in  curvilinear  coordinates.  To  be  more 
precise,  we  put 


US 

Q(V)~ 


p(v)a 

M-s 


9 /i,s 


q(v) 


r3P(v)'1 

(J. 

i8q(v); 


‘(v) 


s>0 


t(V)  s>0 


(36) 


(37) 


and  integrate  the  equations 


.  P* 

Q(V)=J 


(v)  ,  „(v)v^H’s9g 


'(v) 


,  s>0 


(38) 


.(v) 

p 

(IS 


=  I 

M-VM-' 


t3q(v>. 


m 


2P(V)P(V!)  — 


HU’ 


_ +  p(v)p(v)lQ^9V^ 

9(/ 


— -V  — 
.’Q(V)  9Q^  9Q^ 


s>0  (39) 
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together  with  the  equations  of  motion 


P<v)  _ F(V)  .  X  3Vv 


(40) 


aq 


4.  Self-Consistent  Eikonal  Method  (SCENT) 

Let  us  now  introduce  the  average  flow  with  the  kinetic  energy  given  by  the  relation 


with 


1P2«2)-£ZPV<Q) 


P(Q)  =  Ipv(Q)  . 

v 


-*v) 

P  (Q) 


/  p(Q) 


(41) 


(42) 


Multiplying  both  sides  of  Eq.(15)  by  Pv  .  summing  over  v  and  taking  into  account  that  the 
— ^ 

vector  coefficients  F  ( in  Eq.  (17)  change  their  sign  under  the  interchange  of  v  and  v '  one  can 
easily  verify  that  the  function  P(Q)  satisfies  the  equation 

y  P2(Q)  +  V(Q)  =  E  (43) 

with 


V(Q)  =  IVv(Q)pv(Q)/p(Q) 


(44) 


and  hence  we  can  define  the  velocity  field  of  the  average  flow  by  means  of  the  relation 


P(Q)=V  W  , 


(45) 


where  W(Q)  is  Hamilton’s  characteristic  function  for  the  potential  V(Q) ,  namely. 


V  W 


+  V(Q)  =  E  . 


(46) 
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form: 


Let  us  now,  in  following  Micha'  work10,  write  functions  Eq.  (10)  in  the  common  eikonal 


Vv(Q)  =  Xv(Q)exp[iW(Q)/h]  , 


^v(Q)  =  %v(Q)exp(iW(Q)/h], 


where  we  use  the  tilda  to  mark  the  diabatic  wavefunctions  to  distinguish  them  from  those  in  the 
adiabatic  representation.  Let  us  now  neglect  the  difference  between  the  Jacobians  in  the  sum  in  the 
right  side  of  Eq.  (28),  express  the  function  ¥  V(Q)  in  the  right  side  of  Eq.  (29)  in  terms  of 

ft  fr\  T\  — 


^v(Q(v)_t(v),w0.)- Jv  *v  , 


and  substitute  the  resultant  expression  in  Eq.  (28),  writing  the  semiclassical  equations  as 
dp  _  i— >(v)2  — >  ( _»  — »(v')^ 

ih_vL=  v  v-|k  pv+  i  f  •  I. p  +  k  ;p 

dt^.x  L  v  2  Jrv  v,^.v  W’  v 


where  we  put 


-Kv)  -Kv) 

k  (Q)=  P  (Q)-P(Q) 


A  similar  equation  in  the  diabatic  representation  takes  the  foim 


dPv  _  i  ~>(v)2 

sr-‘  K*  ?  h  p  , 

^(v)  V*V  vv  v 


where  H  ,  are  some  diabatic  potential  couplings  and  the  tilda  is  again  used  to  distinguish 
between  the  two  representations. 

The  main  idea  of  the  method  is  to  run  trajectories  on  the  single  potential  surface 
approximating  the  potential  V(Q)  by  means  of  the  relations: 


V(Q)»  !Vv(Q)|Pv(Q)|  /p(Q)  , 
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(54) 


V(Q)«  I  H  J  (Q)P*v(Q)/p(Q)  , 

v  v'  vv  v 

with 


P(Q)  =  S 

v 


Pv(Q 


(55) 


— Kv)  — > 

After  neglecting  the  deviations  k  of  the  momenta  P  from  the  average  P  which  cannot 

be  evaluated  on  the  trajectories  in  question  we  thus  come  to  the  following  equations  for  the 


functions  PV(Q)  and  |5V(Q): 


ift- 


dP, 


dt 


=  [Vv-V]Pv+  J  ..Fw'*ppv" 


V*V 


itIdr1=[Hw-v]Pv+  ?  hvv,pv.. 


v  *v 


(56) 


(57) 


Note  that  both  sets  of  equations  conserve  the  right  side  of  Eq.  (55)  along  any  trajectory  governed 
by  the  potential  V(Q) .  The  equations  are  solved  with  the  boundary  conditions 

P  ^(XWq)  =  P  V«X  wQ)  =  5vl[p°(w0)]1/2exp  (-  ik  jR0  +  iW0(W())/h)  (58) 

Compared  with  the  discussion  presented  in  Section  3,  a  certain  simplification  of  the  formalism 
comes  from  the  fact  that  Eq.  (44)  does  not  contain  the  ratio  of  functions  Ay(Q)  and  Ay,(Q)  and 

hence  one  does  not  have  to  deal  with  the  indeterminate  forms  necessary  to  evaluate  the  effective 
potentials  Vv  at  t=0.  It  is  worth  pointing  out  that  the  neglect  of  the  deviations  of  the  momenta 

— >  ( v)  — ► 

P  from  the  average  P  in  Eqs.  (56)  and  (57)  is  the  most  serious  assumption  made  so  far 
because  it  changes  the  asymptotic  behavior  of  wave  functions  after  the  couplings  between  the 
channels  completely  turn  off.  In  the  next  section  we  show  how  one  can  eliminate  this  defect  by 
running  a  separate  trajectory  in  the  field  of  each  adiabatic  potential,  instead  of  a  single  trajectory  in 
the  field  of  the  average  potential  V(  Q) . 

The  set  of  equations,  Eq.  (57),  has  been  recently  derived  by  Stodden  and  Micha17  by 
evaluating  directly  the  time  derivative  of  the  Jacobian,  Eq.  (24).  It  should  be  emphasized  that 
although  similar  sets  of  coupled  differential  equations  have  been  derived  in  the  literature,  starting 
from  the  time  dependent  Schroedinger  equation  in  either  adiabatic18  or  diabatic19'20 
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representations,  in  the  latter  case  one  propagates  some  time-dependent  coefficients  in  a  basis-set 
expansion  whereas  we  try  to  build  the  wavefunction  itself.  The  crucial  difference  comes  from 

turning  points  where  the  phases  of  the  functions  (5  V(Q)  and  {5  V(Q)  take  a  discontinuous  change. 

When  the  eikonal  method  is  applicable,  it  shoulf  give  more  accurate  results  since  it  treats 
translations  semiclassicallly,  rather  than  classically  as  in  the  time-dependent  expansions. 

The  novelty  of  our  approach  when  extended  to  the  multi-dimensional  case  comes  from  the 
explicit  use  of  the  time-scale  separation  to  carry  those  functions  through  the  caustics.  We  show 
below  that  the  equations  of  motion  used  by  us  in  the  multi-dimensional  case  differ  from  those 
derived  in  the  literature17'19  . 

Potential  Eq.  (54)  has  a  very  interesting  feature,  namely,  its  first  derivative  with  respect  to 
time  is  given  by  the  simple  relation: 

V=  I  P  .  V  H  i3vP*  /p  .  (59) 

v,v 

To  prove  it  we  need  just  represent  Eq.  (54)  and  Eq.  (57)  as 
f* 

Vs  j$  Hf5/p  (60) 


and 

ihP  =  (H-Vl)P  (61) 

and  then  substitute  Eq.  (61)  in  the  derivative  of  Eq.  (60)  with  respect  to  t,  taking  into  account  that 
the  matrix  H  -  V  J.  commutes  with  H  .  The  direct  consequence  of  the  proved  result  is  that  the 
governing  force  in  the  one-dimensional  case  is  given  by  the  simple  expression 

.1  h  dU 

P  =  P  ~ “ f  P  /  P  .  (62) 

-  dQ  - 

(We  can  always  make  g*  *  equal  to  a  constant  l/|i^  by  the  appropriate  change  of  variables  so  that 

no  additional  term  appears  in  the  equation  for  the  streamline  motion.) 

The  governing  potential  force  in  the  multi-dimensional  case  has  a  more  complicated  form: 
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Q  H  F  jj,  =  P 


*■  9  qs 

0/P+  I  q 


9Q 


|i  _ 


s>0 


3Q  L 


2  Re 


(  F 

P  HP, 


K  “H,S 

V  ~  J 


^  h  2 
/P-P.S?  HP/pZ 


(63) 


where  wc  put 


ap 


_  ,  s 


3qS(v) 


ap 


s  9qs 


(v> 


(64) 

(65) 


Evaluation  of  derivatives  Eq.  (64)  along  trajectories  governed  by  the  Ehrenfest  potential  Eq.  (53)  is 
discussed  in  Appendix  B. 


5.  Adiabatic  Velocity  Field  Method  (AVFM) 

As  outlined  in  the  end  of  Section  II  we  propose  to  neglect  quantum-mechanical  correction 
Eq.  (17)  and  to  run  trajectories  on  the  adiabatic  potential  surfaces.  However  to  do  it  one  needs  the 
initial  conditions  for  the  momenta  in  the  channels  which  arc  not  populated  at  the  initial  moment 
Since  kj  is  the  largest  momentum  in  the  system,  these  momenta  can  be  defined  by  the  expression 

hkv=  A r[Vi(R0’wo)  -  Vv(R<Jwo)]  +  h2kl  l66) 


where  the  adiabatic  potentials  represent  cuts  through  the  global  potential  surface  with  fixed 
quasiperiodic  coordinates  establishing  the  asymptotic  limit  For  atom-diatom  collisions  the  different 
channels  correspond  to  excitation  of  the  atom  for  fixed  internal  vibrational  coordinate.  They  may 
also  correspond  to  different  spin-states  of  the  diatomic  molecule,  as  in  the  case  of  NO  scattered  by 
a  silver  surface. 

— > 

We  also  assume  that  dependence  of  the  nonadiabatic  couplings  F  on  the  quasi-periodic 

w' 

12  0 
coordinates  Q  ,  Q  ,...  is  negligible  at  large  values  of  the  streamline  coordinate  R=Q  and  only  the 

component  F^  ,  along  the  streamline  motion  is  asymptotically  important.  As  a  result  of  this 
assumption  SME's.(13)  become  asymptotically  separable.  This  implies  that 
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div  [v  v  p  ]r=oo  =  "  ^vIr  (v  v) 


and  hence 


arg[¥v(°>wo^]  =  _  ikvRQ+  iW  (wQ)/h  +  mod(rc) 

p(}I)(0’wo)=p(^)(0’wo)  H>0  • 


Let  us  now  represent  Eqs.(28)  as 


f^v  ^  r  (  \  i  *4(v-)  2" 

ih  R(V  -)[R’W]  [jr- J W  =  LV  _)[R’w] )  -  E  -  T  p  [  ^  J  QR’W1 


.1/2,  D  . 

J  [R,w] 

-ih  Z  — rjx -  F 

vVvJv?[Q(v-)]  VV 


(  \  i5W  IQ,  .1  /h 

(VjtR.wlJ.  p  [Q(v  )]e  ”L  (V->J  t~I8.w] 


where  minus  super  (sub)-scripts  denote  the  incoming  wave.  The  most  essential  obstacle  in  solving 
these  equations  in  the  case  of  more-than-oue  channels  is  that  the  functions  £~[R,w]  with  v'  ^  v 

in  the  right-hand  side  of  Eq.  (70)  are  supposed  to  be  evaluated  on  the  trajectory  Q(V)[ wq1  • 
(v) 

The  components  [Q(v  of  the  momentum  P  Jacobian  JV'_[Q(V  )]  and 

the  phase  shift  8W  v[Q^y  ^  are  approximated  respectively  as 


f  Ou  (v) 

*  /2P0  P(vj]-P(,} 


\\ F-v/n  \  I  y  Jl‘'p(v'p(v,l2.„2 

[  (Q(v->) '  2  ^0s  P  V  J  P<W 


•(71) 
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(v)  (v) 

p;  [Q(v.)]-p;  [q(v-)] 


1  Ll-1  U  (1+1 

E-V|  R.q'-.Q  ,Q(v),Q  ... 


with 


p  =i  j  g<Vv,) 

P(v1  1^*  H  ■ 

u.'  if  a'  a' 

Q  [R,w]  =  yjQfv  _.[  R,  w]  +  Q(v._}t  R,  w]  J  . 


(vrQ 

’■■■) 

ti>0 

(72) 

(73) 

(74) 

Iv'IQ(v-)]  “  1  v-[Q(v'-)]  „5 


P(vj] 


H20  (V) 

VWv.)] 


(75) 
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(V) 


5W 


(v) 


[Q(V-)]% *>J  ^  8  Wv-)]  • 


(V) 


(76) 


(Sign  (+)  or  (-)  in  Eq.  (72)  is  chosen  depending  on  the  direction  of  the  appropriate  one¬ 
dimensional  motion.)  We  integrate  Eq.  (70)  by  running  a  set  of  trajectories  on  the  adiabatic 
potential  surfaces  with  the  initial  conditions  given  by  Eqs.(68)  and  (69)  for  the  momenta  and 

Cv(a  wQ)  =  5vl[p°(w())]1/2exp  (-  ik  jR0+  iW°(woyh)  (77) 

for  the  wavefunctions  themselves. 

The  most  important  consequence  of  the  presented  mathematical  arguments  is  that  the 
semiclassical  multichannel  solution  sought  for  does  have  the  eikonal  form  --  remember  that 
existence  theorems  have  been  proved3  only  for  single-channel  processes.  The  equations  for  the 
outgoing  wave  are  similar  to  Eq.  (70)  except  that  they  are  solved  with  different  initial  conditions 
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formulated  using  the  WKB  connection  formula  at  the  classical  turning  point  for  the  streamline 
motion  on  the  appropriate  adiabatic  potential  surface. 

One  of  the  most  important  advantages  of  the  new  theory,  compared  with  SCEM,  is  that  the 
wavefunctions  are  propagated  by  means  of  the  single-channel  semiclassical  equations  in  regions 
where  the  nonadiabatic  couplings  vanish.  In  particular,  this  automatically  excludes  population  of 
energetically  forbidden  states,  which  is  usually  observed  if  wavefunctions  are  propagated 
according  to  the  prescriptions  of  SCEM. 

It  is  worth  emphasizing  that  the  most  serious  approximations  are  made  by  us  to  extrapolate 
wavefunctions  from  one  trajectory  to  another  in  the,  space  of  quasi-periodic  degrees  offreedom.  To 
find  a  solution  of  the  one-dimensional  problem  one  just  needs  to  propagate  SME's  (13)  along  the 

~KV)  , 

streamline  coordinate  R  with  P  (R)  evaluated  from  the  classical  trajectory  run  on  the  Vth 
adiabatic  potential  surface,  instead  of  using  the  exact  self-consistent  potential  Eq.  (16). 

6.  Evaluation  of  the  S-matrix  elements  from  semiclassical  wavefunctions 

To  avoid  the  singularities  at  the  classical  turning  points,  when  evaluating  the  S-matrix 
elements: 

( v) 

(v)  —  ik  ,  R. 

k  e  v  I  Jdw\j/W(R  w,u)u  (w)  .  (78) 

v  n=+i  1 


where  uy,(w)  is  an  exact  quantum-mechanical  wave  function  and  v  is  a  set  of  the  appropriate 
quantum  numbers,  we  first  represent  Eq.  (78)  as 


,(v) 


Rf+AR 

j  dR 
v;  v  AR  DJ 

Kf 


.(v) 

*  v;v' 


(79) 


The  relation  turns  into  the  identity  if  we  deal  with  the  exact  matrix  element  which  is  independent  of 
R.  We  thus  find 


AR 


(v)_(v) 
v  v;v’ 


I 

u=±l 


t^w^u) 


/Tv)  f  j/2  s 

Idxjk  (wQ,u)  }  dt  Jy  (t,w0,o)C  vv(t,w0,-u)uv,(w(t;w()))  , 


(80) 


31 


i  2 

where  dx  =  dx  dx  ... , 


C w(t,w0,1))  ~  ^w(t,wo,^)< 


(V) 


-ik  (w0,-o)R(t,w0,\)) 


(81) 


=  ps(ws)  /m  >0 

dx5*  s  (82) 


(with  ps  and  ms  used  for  the  momentum  and  the  effective  mass  of  the  sth  quasi-periodic  degree  of 
freedom)  and  Ty  is  the  absolute  value  of  determinant  of  the  auxiliary  matrix  with  the  elements 


P(v)(wff')> 


(83) 


f-ts 

propagated  instead  of  the  functions  to  avoid  singularities  in  the  initial  conditions  for 

functions  Eq.  (37).  (The  function  P^(wq)  in  Eq.  (77)  is  chosen  in  such  a  way  that  the  matrix 


element  is  equal  to  1  for  elastic  collisions.)  The  crucial  advantage  of  the  integral  Eq.  (80) 
compared  with  Eq.  (78)  is  that  the  Jacobian  J  v  appears  now  in  the  numerator  and  it  is  multiplied 
by  function  Eq.  (8 1)  which  is  nonsingular,  in  contrast  to  the  function  \|/  w  (see  Eq.(29) ). 


7u-C9nglu?iQp§ 

The  analysis  presented  in  this  subsection  revealed  some  promising  new  directions  in 
practical  implementation  of  multichannel  multidimensional  semiclassical  theory.  We  have 
developed  a  self-consistent  approach  that  makes  it  possible  to  address  (under  some  reasonable 
assumptions)  the  most  important  issues  of  semiclassical  theory.  The  following  computationally 
affordable  algorithms  could  be  cited  as  the  main  results  of  this  subsection: 

(1)  well-defined  prescriptions  for  continuation  of  global  (multidimensional)  semiclassical 
wavefunctions  through  classical  turning  points, 

(2)  a  practical  numerical  algorithm  to  extract  final  populations  of  vibrational  levels  from  the 
semiclassical  outgoing  solution, 

(3)  implementation  of  the  single-surface  developments  (1)  and  (2)  in  the  multichannel 
theory  with  valid  use  of  the  Ehrenfest  effective  potential  (SCEM),  and 

(4)  development  of  the  new  multichannel  semiclassical  theory  referred  to  as  "Self- 
Consistent  Adiabatically-Corrected  Eikonal  Method"  (SCACEM)  that  goes  beyond  the 
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assumptions  underlying  SCEM  and  explicitly  takes  into  account  specific  topology  of  the 
adiabatic  potential  surfaces. 

We  have  obtained  some  preliminary  but  very  encouraging  numerical  results21  in  support  of 
our  treatment  of  the  single-surface  problem  by  considering  vibrational  energy  transfer  in  atom- 
diatom  collisions  in  the  IOS  approximation22.  Another  project  which  was  recently  initiated  to  test 
the  validity  of  the  developed  formalism  is  a  study  of  rotational  excitation  of  the  rigid  NO  molecule 
scattered  by  the  uncorrugated  silver  surface.  The  main  atraction  of  this  problem  is  that  exact 
quantum-mechanical  calculations  have  been  performed  by  Smedley,  Corey  and  Alexander23, 
including  electronically  inelastic  scattering  due  to  nonadiabatic  couplings  between  different  spin- 
orbit  and  A-doublet  states.  Excellent  agreement  between  those  calculations  and  the  one¬ 
dimensional  reduction  of  SCEM,  treating  rotations  quantum-mechanically  has  been  reported  in  an 
earlier  publication24;  our  next  step  is  to  construct  semiclassical  wavefunctions  from  realistic 
classical  trajectories:  this  is  being  done  for  NO  approaching  the  smooth  rigid  surface  to  compare 
the  obtained  state-to- state  information  with  the  results  of  the  previous  studies. 


Appendix  A 


Let  us  consider  nutational  motion 
Schroedinger  equation: 


of  a  diatomic  molecule  described  by  the 


2  1  d 
sin  0  d0 


mV 

sin  ^0 


2I*(V(e>'ew)}'n.  =  0 


(84) 


where  the  angle  0  is  restricted  to  the  interval  (0,7t).  We  consider  the  general  case,  allowing  the 
molec  jle  to  rotate  around  the  quantization  axis.  We  will  need  such  an  extension  when  considering 
the  rotationally  inelastic  scattering  of  diatomics  by  a  corrugated  rigid  surface,  because  molecule- 
surface  interactions  couple  states  with  different  m.  In  following  Langer’s  suggestion15,  we  should 
now  transform  the  finite  interval  (0,7t)  onto  the  whole  axis  (  -  oo,+  «>) .  As  initially  emphasized  by 
Froman  and  Froman25  such  a  transformation  does  not  lead  to  a  unique  modification  of  the  Bohr- 
Sommerfeld  quantization  rule.  To  eliminate  the  ambiguity,  Adams  and  Miller26  suggested  choosing 
the  transformation  that  reproduces  the  exact  quantum  spectrum  of  free  motion.  This  suggestion 
seems  especially  appropriate  for  the  scattering  problem  in  question.  We  show  below  that  the 
mapping  transformation 

x=ln(tg|).  •  (85) 
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leads  one  to  the  semiclassical  energy  levels  shifted  up  by  the  same  constant  h^/(8I*)  from  the 

exact  rotational  levels  j(j+l)  h^/(2I*) .  The  appropriate  semiclassical  wavefunctions  have  been 
discussed  in  detail  by  Landau  and  Lifshitz  who,  however,  come  to  this  result  only  after  neglecting 
the  singular  quantum-mechanical  correction,  whereas  we  prove  that  transformation  Eq.  (85) 
mapping  the  interval  (0,;t)  on  the  whole  axis  (  -  <»,+  oo)  simply  eliminates  it.  In  fact,  taking  into 

account  that 

ex  =  tgy  (86) 

we  find 

dx  _  1 

d9  “  sin  0  (87) 


and  hence  Eq.(84)  takes  the  form: 


■^^r4(mh)2+v(x) 

L  2  dx2  2 


¥  =0 
m 


(88) 


with 


v(x)  —  -  4I*[ V(9(x))  -  e]  e2x  (l  +  e2x) 


-2  . 


(89) 


We  thus  find  that  -  (mh )  /2  now  plays  a  role  of  the  energy  levels  in  the  potential  well  v(x).  In  the 

2  iff 

particular  case  of  free  nutational  motion  V(  9 )  =  0  and  -  I  £  £  v(x)  <  0,  i.  e.,  -  (mh )  >  -  21  e 
and  hence 


e> 


m2h2 


(90) 


Applying  the  usual  WKB  rules  to  potential  (89): 


j 

j  dx-J -  (mh)^  -2v(x)  =  7th(n  +■  1/2) 
x  _ 


(91) 
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and  coming  then  back  to  the  angle  8  we  find  that  the  semiclassical  energy  spectrum  for  Eq.(84) 
is  given  by  the  relation: 


w-v®)- 


mV 

sin  ^0 


=  7*(n.+  1/2)  . 


(92) 


For  free  nutational  motion  the  change  of  variables  r=cos  0  leads  one  to  the  integral  easily 
evaluated  by  means  of  the  method  of  residues,  by  analogy  with  a  similar  derivation  given  by 
Sommerfeld  for  the  Coulomb  problem  (see  Ch.9  in  Ref.  12  ).  The  residues  at  the  points  -1,  +1  and 


infinity  are  equal  to  R  ^=ilml/2,  R+^=ilml/2  and 


R~=- 


-"'J 


2I*e. 


jjmi 


respectively  so  that 


-Imlh 


=  h(n  +  1/2) 


(93) 


and  hence 

_  h2j(j  +  l)  h2 

£j,lmr  21*  +  81* 


(94) 


We  see  that  quantization  condition  Eq.  (92)  reproduces  the  exact  values  for  rotational  frequencies 
but  the  zeroth  point  energy  turns  out  to  be  larger  by  h2/(8I*)  . 


Appendix  B. 


Evaluation  of  the  forces  according  to  SCEM  recipe  Eq.  (63)  is  complicated  by  the  fact 
that  the  set  of  equations 
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0Q_ 

d  ws 


=  1 


tni  H  sag 

z  P  +  P  10  — - 
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r  C 

is  unclosed  because  we  need  to  know  the  derivatives  of  the  functions  Q  with  respect  to  q  : 


M-s  2  4 

3Q  _  3  Q  4 

aqs'  aqsaqs'  ss' 


(100) 


as  follows  from  the  relation: 


3  f  acl  1-  v  a<l  ^  3cl 
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(101) 
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obtained  by  differentiating  the  identity 
8qs  \iY 

I— =8  ,  „ 

,  u.  ss 

M-  dQ 


(102) 


with  respect  to  qs.  The  situation  is  completely  different  from  that  for  the  single-channel  problem 
because  in  the  latter  case  the  potential  is  a  known  function  of  the  coordinates. 

|i' 

To  overcome  this  difficulty  we  evaluate  the  functions  Qss>  in  Eq.  (100),  assuming  that  the 
problem  is  adiabatically  separable,  namely,  differentiating  the  equation: 


ws(t,w*) 

or  S  c  S 

'  1  dw  /  ps(w  ) 

ws(trwS) 


t  -  t^  =  m 


(103) 


with  respect  to  w^  we  find 


<r,(t,ws0)/Ps(t,ws0)  *  Q^(t1,ws0)/Ps(t1,w^)) 


(104) 


and  hence 


.  tis 

Q  *VV/FV 


(105) 


Differentiating  Eq.  (102)  with  respect  to  w^  we  conclude  that  propagation  of  the  functions 

M- 

Q  .  in  time  is  governed  in  the  aforementioned  approximation  by  the  equations: 


C=^s[v'vvv/p;;]. 


(106) 
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He-He*  Scattering  Calculations 


L  Model  Validation  Calculations 

Quantal  Scattering  Validation  of  the  SCE  Methodology 

The  validity  of  the  SCEM  for  describing  the  quenching  of  electronically  excited  atoms  in 
the  gas  phase  has  been  tested  by  comparison  with  quantum  mechanical  calculations  on  the  He2 
system.  This  model  is  based  on  an  analysis  of  the  relevant  electronic  states  of  the  He2  system  in 
the  Q  basis. 
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Coupling  (eV) 


We  examined  a  model  which  includes  only  coupling  by  nonadiabatic  processes:  the 
quenching  of  the  2p*P  state  of  He  to  the  2s3  S  state  in  collisions  with  a  ground  state  He  atom.  This 
model  includes  the  2  -  1  g  ,2-0"  and  3  -  1  g  molecular  states.  The  2  -  1  g  state  correlates 
asymptotically  with  the  2p*P  atomic  state  and  has  the  B1^  AS  state  as  the  dominant  component. 
The  2  -  0g  and  3  -  1  g  states  are  nearly  degenerate,  correlate  asymptotically  with  the  2s3S 
atomic  state,  and  have  the  c3Ig+  AS  state  as  the  dominant  component.  Denoting  these  molecular 
states  1,  2,  and  3,  respectively,  the  electronic  Hamiltonian  in  a  diabatic  Q  basis  is  given  by 


H  =1 

el 
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2(J.R 

SO 


13 
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2|lR 

V. 
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2p.R" 


SO 


13 


23 


2|iR 

V, 


007) 


where  V1  and  V2~V3  are  the  adiabatic  potentials  for  the  B  and  approximately  degenerate  c  states, 
respectively,  p.  is  the  reduced  mass  of  He2,  J12  and  J23  are  matrix  elements  of  the  electronic 
angular  momentum  ladder  operator  between  the  state  pairs  1-2  and  2-3,  respectively,  and  S013  is 
the  matrix  element  of  the  spin-orbit  operator  between  states  1  and  3.  For  the  present  calculations 
the  adiabatic  potentials  are  approximated  from  known  spectroscopic  information.  The  couplings 
are  approximated  as  Gaussians  centered  around  the  crossing  point  between  the  B  and  c  states. 

The  relevant  experimental  features  to  be  calculated  typically  include  the  shape  of  the  cross 
sections  for  electronic  transitions  as  a  function  of  the  collision  energy.  Figure  2  shows  a  favorable 
comparison  between  the  SCEM  final  cross  section  results  for  this  model  with  exact  quantal 
scattering  results.  Keeping  in  mind  that  the  SCEM  is  a  semiclassical  theory  of  electronic 
transitions,  the  excellent  agreement  of  Fig.  3  for  detailed  opacities  (solid  line:  EXACT;  dashed  line: 
SCEM)  provides  a  positive  step  in  the  validation  of  the  SCEM  technology  being  employed  in  this 
research. 

Being  a  semiclassical  description,  the  SCEM  framework  is  typically  employed  to  describe 
either  a  forward  process  or  the  reverse  process  of  a  quantum  mechanical  state-to-state  transition  but 
not  both  within  the  same  computation.  Achievement  of  state-to-state  lime  reversal  can  however  be 
readily  checked  by  making  comparisons  of  the  SCEM  treatment  of  excitation  and  quenching 
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0.02  0.03  0.04 


Translational  Energy  Etr  (au) 

Figure  2.  Cross  section  for  He^P)  quenching  to  He(3S).  The  comparison  is  between  exact 
quantum  results  and  semiclassical  SCEM  results  within  the  infinite-order- sudden 
approximation  (IOS). 


Impact  Parameter  b  (a.u.) 

Figure  3.  Opacity  function  for  He^P)  quenching  to  H e(3S).  Methods  are  as  in  Fig.  2. 
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dynamics  separately.  This  has  been  done  for  the  He2*  system  by  employing  detailed  balance  as 
follows: 


,2  .2 

k.a.  .  =  k.a. 

‘  »  -»J  J  j-*» 


(108) 


where 
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1  ^  tr,  i 
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h  h 

Detailed  balance  was  tested  for 

was  found  that 


-  W.) 

»'  (109) 

Etot=0.82833H,  W1=0.72833H,  W2=0.77972H  and  it 


k2lal_^2  =  l  127.7  and  k ]a2  ^  (  =  1238.5 

thereby  displaying  an  excellent  10%  agreement  of  detailed  balance  from  these  completely  distinct 
SCEM  runs. 


2-  Eadiativg.,vgrs»s,.Nonradiativg  Quenching  in  the  gas  Phasg 

Theoretical  Treatment  of  the  Photon  in  Radiative  Quenching: 

We  treat  the  problem  of  radiative  quenching  in  the  same  formal  framework  as  used  on  an 
earlier  occasion  to  treat  photoexcitation  in  polyatomics.  This  is  feasible  in  principle  because  these 
processes  are  inversely  related  to  each  other  and  involve  the  same  underlying  hamiltonian;  also, 
this  is  conceptually  useful  because  of  the  structure  of  the  SCE  methodology.  We  only  briefly  point 
out  certain  essential  aspects  of  this  framework  to  the  problems  now  on  hand  and  the  validation 
requirements  of  the  SCEM. 

The  radiative  quenching  hamiltonian  within  the  electron-field  representation  can  be  written 
as: 


H  =  -  (h2  /  2\i)  V2  +  H  .  +  H  +  H. 

^  R  el  rad  int 


with 


(110) 
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H  .  =  H  .  +  H 
el  el  s.  a , 

.(0) 


(111) 


where  H  d  contains  the  electronic  kinetic  energy  as  well  as  the  interactions  of  the  electrons  and 
the  nuclei  among  themselves  whereas  H  Q  denotes  the  spin-orbit  interaction  which  often  plays 
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an  important  role  in  the  stability  and  collision  dynamics  of  HEDM's.  We  will  employ  a  second 
quantized  form  of  the  framework  with 

H  j  =  X  CaCaEa(R) 

“  (112) 

denoting  the  diagonalized  operator  He[i  and  Ea(R)  the  energies  of  the  kets  |  a  )  |  of ).  .  . 

representing  the  resulting  adiabatic  electronic  states.  Ca  and  C*  therefore  denote  state  annihilation 
and  creation  operators.  The  other  terms  in  the  hamiltonian  are  the  multidimensional  Laplacian  for 
the  collisional  and  internal  nuclear  state  (projectile  and  target)  nuclear  kinetic  energies,  the  radiation 
field  hamiltonian 


H  =  y. h  v,  a+  a. 
rad  r*  ka  ka  kcj 
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and  the  radiation-matter  coupling  term 
^  hv,  ^ 
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written  within  the  electric  dipole  approximation.  It  is  seen  that  is  controlled  by  the  magnitude 
of  the  electronic  matrix  elements  of  the  molecular  total  dipole  operator 


_eDaa<R)  =  > 


(115) 


with  tid  —  eXr:  +Z  Z,  R , ,  where  r,  denotes  all  electron  coordinates,  RT  all  nuclear  ones, 
i  I 

e  is  the  absolute  value  of  one  electronic  charge  and  Zj  represents  the  nuclear  charge.  ak(j  and 
a£ff  are  photon  annihilation  and  creation  operators  respectively  for  second  quantized  radiation 

field  modes  in  a  volume  V.  e0  denotes  the  vacuum  electric  permittivity  and  is  a  unit  vector  of 

the  polarization  of  the  radiation  (electric)  field.  With  these  definitions,  the  time-independent 
Schrodinger  equation  for  the  radiative  quenching  from  an  electronically  excited  polyatomic 
molecule  that  is  isolated  or  in  a  collisional  process  is  given  by, 

HI¥(R)>  =  E(R)  ly(R)  >  , 

e,f  e,f  (H6) 
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E(R)  denoting  the  total  energy  of  the  electron-field  wavefunction,  I  \j/(R) )  for  nuclear 

C,  I 

configuration  R,  which  is  in  turn  readily  expanded  as 


I  V(R)>e 
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in  terms  of  products  of  the  adiabatic  spin-orbit  states  j  a  )  with  photon  number  states  |  N  ) 

ko 


yielding  the  infinite  set  of  coupled  equations  for  the  nuclear  wavefunctions  n  (R); 
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H 

where  the  hamiltonian  matrix  qn  a  •  n'  includes  the  various  intra  and  inter-  electronic  state 

ka  k'a’ 

dipole  couplings  and  momentum  (i.e.,  radial  derivative  and  rotational)  couplings.  The  problems  to 
be  treated  below  are  now  readily  perceived  as  simple  cases  of  this  general  radiative  quenching 
hamiltonian  that  are  pertinent  for  the  subtask.  A  finite  set  of  equations  will  arise  from  truncating  the 
hamiltonian  matrix  to  include  just  the  relevant  channels  for  the  computations. 


Calculations  with  accurate  potentials  and  couplings: 

Based  on  the  recent  calculations  of  accurate  potentials  and  couplings  by  Yarkony,1  we  have 
performed  several  dynamics  calculations  which  are  described  below. 

We  examine  the  pathways  of  collisional  quenching  of  He(23S)  which  will  become  more 
important  at  higher  gas  phase  pressures.  We  assume  only  binary  collisions  are  important  and 
quenching  occurs  by  nonadiabatic  dynamics  involving  scattering  on  multiple  electronic  surfaces, 
notably  the  a(3Iu+),  A(1IU+),  b(3ng),  B^TIg),  c(3Zg+),  and  C^Ig4-)  states  of  He2. 

Electronically  nonadiabatic  collisions  of  He^S)  and  He(23S)  atoms  can  begin  on  either  the  a  or  c 
state  potentials  and  ultimately  lead  to  quenching  to  the  X-state.  The  a(3Lu+)  and  c(3Zg+)  states  of 
He2  correlate  asymptotically  with  He('S)  (ground  state)  +  He(23S).  Dynamical  studies  of  the 
important  pathways  leading  to  such  an  overall  spin-forbidden  quenching  can  now  be  studied 
employing  Yarkony’s  new  ab  initio  results1  for  the  relevant  potential  energy  surfaces  and 
nonadiabatic  couplings  of  this  metastable  system.  Since  these  participating  states  involve  a  mixture 
of  singlets  and  triplets,  transitions  are  to  be  described  in  the  total  angular  momentum  (£2)- 
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Table  2:  He2  electronic  states  in  the  &  basis 


Q  state 

degeneracy 

spectroscopic 

notation 

in  A-S  basis 

dominant 

space-spin 
configuration  (Ms) 

asymptotic 
atomic  states 

i-o; 

1 

X 

‘r,(0) 

ls2lS,  ls2lS 

1-0; 

1 

a 

3r»  (0) 

ls2s  3S,  Is2  'S 

1-1. 

2 

a 

3ra(i),3i+tt(-n 

ls2s  3S,  Is2  'S 

i-o; 

1 

A 

'z;  (0) 

ls2s  ’S,  Is2  ‘S 

2-0; 

1 

b 

3nI.1(-i)-3n1.-i(i) 

ls2p  3P,  Is2  'S 

i  ~  o  g 

1 

b 

3ngl(-i)+3ng_1(i) 

ls2p  3P,  Is2  *S 

1-1* 

2 

b 

3nfcl  (0),  3ng  (0) 

ls2p  3P,  Is2  *S 

1-2, 

2 

b 

3ng_1(-i),3ngl(i> 

ls2p  3P,  Is2  'S 

2-1* 

2 

B 

lng  .,(0),  ‘ng  t(0) 

ls2p  'P,  Is2  'S 

2  -  0- 

1 

c 

3z;  (0) 

ls2s  3S,  Is2  'S 

3-1. 

2 

c 

3z;  (i),  3r,(- 1) 

ls2s  3S,  Is2  ’S 

3  -  0, 

1 

C 

lK  (0) 

ls2s  'S,  Is2  :S 

Coupling  selection  rules: 
dipole  coupling,  p: 
radial  derivative  coupling,  8/9R: 
angular  momentum  coupling,  J+: 

g«-»u,  A£2  =  0,±1 
g<-»g,  u<->u,  AQ  =  0 
g*-»g,  u«->u,  A£2  =  ±1 
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Table  3:  Total  coupling  in  the  £  basis 


Table  3  (continued)  :  Total  coupling  in  the  ft  basis 
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representation.  Tables  2  and  3  indicate  the  resulting  states  along  with  the  appropriate  radiative  and 
nonadiabatic  couplings. 

Figure  4  illustrates  the  important  triplet-triplet,  singlet-triplet  and  singlet-singlet  dynamical 
pathways  for  the  a  —>  X  quenching  process.  (Please  note  that  the  singlet-triplet  diagram  suppresses 

intramultiplet  couplings).  The  solid  arrows  indicate  transition  dipoles  and  electronically 
nonadiabatic  dynamical  couplings  known  from  Yarkony's  calculations  (Hso:spin-orbit,  (J.sf.spin- 

forbidden  dipole,  and  L:  rotational  couplings)  and  the  dashed  arrow  is  a  negligible  (due  to  large 
energy  gap)  Hso  coupling.  The  diagrams  in  Fig.  4  can  be  employed  to  identify  two  key  pathways 
of  quenching:  (1)  a  direct  spin-forbidden  radiative  pathway  between  the  a-state  and  an  X-state 
dressed  with  a  photon2  (as  seen  in  the  middle  diagram  of  Fig.  4)  and  (2)  a  more  complex 
multistate  dynamical  route  involving  hops  among  a  variety  of  triplet  states  and  eventually  crossing 
into  the  singlet  manifold,  where  strong  quenching  propensity  prevails. 

Although  the  system  can  in  principle  quench  nonradiatively,  it  was  found  that  the  related 
probabilities  are  very  low,  ~10'16,  compared  to  ~10r&  for  the  radiative  routes. 

The  direct  spin-forbidden  radiative  pathway  is  studied  using  a  two-state  model  using  the  a 
and  X  states  within  the  self-consistent  eikonal  method  (SCEM).  The  relevant  potentials  and  the 
spin  forbidden  radiative  coupling  are  shown  in  Fig.  5  where  the  X  state  has  been  shifted  up  by  the 
relevant  atomic  transition  energy.2  Calculations  using  SCEM  and  based  on  Eq.  (118)  yielded  the 
following  cross  sections  (in  aQ2)  for  collision-induced  radiative  quenching:  1.23xl0'4 

(probability- lx  10'6),  1.92xl0‘13  (probability  -l.xlO"7)  at  translational  energies  5.4  and  .27  eV 
respectively.  Probabilities  drop  down  to  1.x lO®  at  .004  eV  but  climb  again  to  l.xlO'7  at  .002  eV. 
These  numbers  are  consistent  with  the  experimental3  cross  section  upper  limit  of  ~l.xlO"9a02 

which  dealt  with  pressure  effects  at  low  translational  energies.  It  is  also  seen  that  there  is  some 
room  for  additional  quenching  contributions  while  remaining  consistent  within  the  experimental 
upper  bound. 

Several  different  models  of  the  multistate  pathways  were  studied.  So  far,  our  results 
indicate  the  following  notable  radiative  quenching  features.  In  spite  of  the  many  additional 
radiative  and  nonradiative  pathways  available  in  the  many-state  model,  the  nonradiative 
nonadiabatic  route  to  the  singlet  manifold  is  incapable  of  competing  with  the  propensities  of  the 
above  two-state  radiative  route.  This  is  consistent  with  Yarkony's  anticipation1  of  a  bottleneck  in 
the  multiplet  crossover  rate.  The  dynamical  problem,  however,  becomes  very  complicated  when 
the  various  photon  emitting  pathways  are  included.  A  full  study  requires  a  thorough  scan  of  photon 
frequency  dependence  in  each  dipole  coupling:  there  are  13  dipole  couplings,  9  allowed  and  4 
forbidden.  Systematic  examination  of  selected  many-state  models  and  trajectories  were  done  to 
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Figure  4.  He2  Transition  Pathways,  based  on  ab  initio  results  of  Ref.  1. 
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Triplet-Triplet  transition  pathways 


Singlet-Triplet  transition  pathways 


Singlet-Singlet  transition  pathways 
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Figure  5.  He2  two- state  (a  and  X)  radiative  quenching  model  potentials. 
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examine  the  sensitivity  to  photon  frequencies  and  dependence  of  results  on  whether  keeping  or 
dropping  the  radiative  couplings.  Althouh  dramatic  changes  have  been  found  in  the  relative 
propensities,  they  remain  relatively  small  compared  to  the  direct  route.  Hence  it  is  probably  safe  to 
a  priori  rule  out  competition  from  these  pathways  at  this  stage. 
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Accurate  Computations  of  Energy  Transfer  in  Atom-Diatom 

CoIlisions:Comparative  Study  of  Vibrational  Excitations  using 

Semiclassical  Wavefunctions  and  Coupled  Quantal  States 

1.  Introduction 

Molecular  collisional  energy  transfer  presents  a  computationally  challenging  problem  to  the 
dynamicist  even  after  the  global  potential  energy  surface  is  known.  A  rich  variety  of  quantum 
mechanical,  semiclassical  and  classical  models  may  be  constructed  to  solve  the  problem  depending 
on  the  behavior  of  the  potential  energy  surface  which  controls  the  collision  dynamics.  The 
information  from  modem  laser-based  experiments,  involving  precise  microscopically  detailed 
observations  of  collisions,  is  in  the  form  of  state-resolved  cross  sections  and  demand  an 
explanation  from  theory  at  that  level  of  detail.  Classical  mechanics  in  itself  does  not  offer  a  long¬ 
term  solution  for  obtaining  such  state-resolved  information.  Fully  quantum  mechanical  treatments 
become  unwieldy  for  realistic  systems  despite  recent  triumphs  that  exploit  efficient  algorithms  on 
supercomputers.  Thus  there  is  a  present-day  need  to  develop  and  exploit  semiclassical  dynamical 
theories  that  hold  promise.  In  order  to  obtain  general  results,  it  has  to  be  based  on  a  systematic 
process  involving  detailed  benchmarks  against  accurate  quantum  mechanical  results  for  smaller 
systems  and  the  development  of  controlled  computational  schemes  which  are  careful  paths  to 
accomplish  the  devious  transcription  from  classical  trajectories  to  semiclassical  wavefunctions. 

The  present  subsection  reports  a  study  of  vibrational  excitation  during  atom-diatom 
collisions  in  gas  phase  and  is  one  of  many  successful  numerical  experiments  in  which  we  have 
employed  the  eikonal  ansatz  for  the  semiclassical  wavefunction.  The  relevant  general  formal 
discussions  have  been  presented  in  the  proceeding  section  in  greater  detail.  The  present  section  will 
focus  on  the  formulation  for  the  vibrational  excitation  problem  and  its  numerical  investigation. 

In  this  subsection  we  consider  a  two  degree  of  freedom  model  of  atom-diatom  collisions 
(O+HF  and  Na+H2  collisions  are  studied)  as  defined  in  the  rotational  infinite  order  sudden  (IOS) 

approximation.  The  numerical  tests  presented  here  are  benchmarked  against  previously  existing 
quantum  IOS  (QIOS)  results.  Two  distinct  dynamical  models  are  employed  in  this  subsection  to 
describe  state-to-state  vibrational  dynamics:  (1)  is  based  on  combining  a  classical  treatment  of 
translational  motion  with  a  quantal  states  expansion  in  target  vibrational  states  (called  RS  VE:  for 
rotational  sudden  with  vibronic  expansion)  and  (2)  is  based  on  employing  semiclassical  eikonal 
wavefunctions  for  all  degrees  of  freedom  obtained  from  their  classical  trajectories  (called  SIOS:  for 
semiclassical  infinite  order  sudden). 
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2.  Theory 

1.  IOS  Approximation  in  the  diabatic  representation 

The  rotational  IOS  approximation1  for  the  atom-diatom  scattering  problem  is  based  on 
solving  the  two-dimensional  Schroedinger  equation  with  the  Hamiltonian: 


A 

h 

IOS 


/ 


9R2 


h2  32  (  h2l(l  +  l) 

2^BC  ^r2  2(lR2 


h2j(j  +  l)" 
2^BCf2  > 


1 +  H(R,r,y) 


(119) 


where  H(R,r,y)  is  the  n  x  n  matrix  of  the  diabatic  potentials  and  couplings  between  different 

electronic  states,  R  is  the  distance  between  the  center  of  mass  of  the  diatomic  BC  and  the  atom  A,  r 
is  the  bond  length  in  the  diatomic  BC  and  7  is  the  angle  between  the  diatomic  and  the  radius 
drown  from  its  center  of  mass  to  the  atom  A. 


2.  Quantum-mechanical  treatment  of  vibrations. 

In  following  Redmon  et.  al.  2  we  reduce  the  problem  to  the  infinite  set  of  the  ordinary 
differential  equations  of  the  second  order  by  using  the  matrix  representation  for  the  hamiltonian  in 
the  space  of  the  vibrational  coordinate  r,  namely,  we  use  the  eigenfunctions  of  the  Schroedinger 
equation: 


h2  d2 


2^1F  +  Vbc«  + 


k2j(j  +  1)  j 
— - $ —  e» 


uJv  =  0 


(120) 


as  the  basis  set  and  represent  Eq.  (120)  as 

dV 


nv 


dR" 


'i^v  +  nlh„V.n'vMVnV  =  0- 

n  ,v 


(121) 


where  we  put 

jl 


x  *  *  Kl+1) 

h  .  ,(R,y)  =8  ,5  . - z — 

nv,nV  nn  w  _2 


+  fl{-  +  ,dr  “vW  -  5„„’vBc(r)]}  • 


(122) 
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with 


*%  ■  '  ev)] 


1/2 


(123) 


and 


AD^  *  Hn(R,r,y)] 


(124) 


used  to  evaluate  the  probabilities  of  transitions  from  the  TJth  electronic  state. 

We  solve  the  set  of  coupled  differential  equations  Eq.  (121)  by  representing  wavefunctions 
Vj,v  in  the  common  eikonal  form3: 


jl  1 

Vnv(^T)  = 


■PnV(R;Y)exp[iWjl(R;Y)/h] 


P  (R;Y) 


(125) 


where  (R;y)  is  Hamilton's  characteristic  function4  found  from  the  Hamilton- Jacobi 

.  4 
equation  : 


2m- 


jl 

P  (R;y) 


2  -Jl  12 

+  V  (R;y)  = 


(126) 


with 


J1  _  dW 

dR 


(127) 


and  V  (R;y)  denotes  the  Ehrenfest  potential: 


J  J 


V  (R;Y)  =  „vJvh-nvPnvPnv/I 


P' 


nv 


(128) 
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The  functions  pnv  satisfy  the  set  of  coupled  equations: 

J  jl  jl  _  jl  jl 

lh3nv  =  ^  .  ,hnv  nV'^nV' ~  V  ^nv 

nv,nV  uv,uv  “ v 


and  are  solved  together  with  the  equation  of  motion  for  the  potential  Eq.  (128). 

3.Semiclassical  treatment  of  vibrations 
A.Equations  of  Motion 

We  represent  the  wavefunction  V(R»  r,  Y4j )  in  the  eikonal  form: 
\|/(R,r,y,lj)  =  A(R,r;y,lj)  exp  [iW(R,r;y;lj)/h] 


where  W(R,r;y;lj )  is  Hamilton's  characteristic  function 

1  ( dW  \2  1  f  dW\2  rT/n  _ 

27TVdrJ  +2ZVdr)  +U(R,r;y;lj)  =  E 


2H 

for  the  potential 


h2l(l  +  1)  h2j(j+l) 
U(R.r,r,lj)  =  V(R,r;y)  + - + 


2nR 


2^BCr2 


(129) 


(130) 


(131) 


(132) 


(133) 


The  momenta  conjugate  to  the  coordinates  R  and  r  are  given  by  the  relation: 

dW  dW 

P0(R,r;y;lj)  =  P1(R,r;r,lj)  = 

To  cany  the  wavefunction  V(R,r,y,lj )  through  turning  points  we  represent  it  as 
-1/2 

V(t,r0;y,lj)sJ  (t,r0;y,lj)  A(0,r0;j)  C(t,r0;r,lj)  (134) 

where  t  is  time,  Tq  is  the  initial  value  of  r  on  the  trajectory  run  on  the  potential  surface  U  and 
J  (t,r^;  ylj )  is  the  Jacobian  of  the  transformation: 


t,rQ  — » R(t,r0;y;lj),  r(t,rQ; y;lj) 


(135) 
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and  then  assume5  that  the  function  ^(t.rgjy.lj)  is  continuous  everywhere  by  analogy  with  the 
one-dimensional  case.  We,  in  following  Ref.5,  solve  the  equation: 


U-E- 


(po)2 

24 


(Pl): 


2  4 


BC 


(136) 


as  if  the  function  ^(t,  r^; y;lj )  itself  were  continous,  finally  shifting  its  phase  by  rc/4  times  the  total 
number  of  the  turning  points  passed. 

B.  Initial  conditons 

Let  us  put  =R  ,  =r,  w*  =  r 


h2j(j  +  1) 

Vj(r>3VBc(r>+V7T 

^BC 

Then  the  initial  conditions  take  the  form 


(137) 


Q°(0,r0,±l)sR0,  Q1(t0,r0,±l)=r0 
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(138) 
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(140) 
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C.  Connection  formulae 

If  the  interaction  between  vibrational  and  collisional  motions  is  negligible  the  derivative  of  r 
with  respect  of  its  initial  value  r^  can  be  evaluated  by  differentiating  the  equation: 


r(t,rQ) 


t  =  nBc  _f*'/pr(f) 


0 


(148) 


with  respect  to  r^  at  the  fixed  moment  t  We  find 

p7l<r{^).'p;1<r°)=0 


(149) 


and  hence 
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V,  0  7 


=pr(rVpr(ro) 
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Therefore  if  the  off-diagonal  element  Q  ^  can  be  also  neglected  the  Jacobian  J  factors  as 
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(151) 


|jJ (t»r 0»±  l)  =  P°(t,rff±  l)P1(t,rff±  D/P1  (tQ,rff  ±  1)  , 

and  its  zeros  approximately  coincide  with  zeros  of  the  momenta  P^  and  P*.  We  thus  assume  that 
jumps  of  the  Jacobian  on  the  caustics  are  governed  for  adiabatically  separable  motions  by  the 
standard  WKB  rules,  namely,  that  the  function  l£(Q)t  is  continuous  there,  whereas  the  phase  of 
the  function  C(Q)  decreases  by  rc/2: 

W(Q)— >W(Q)-«/2  (152) 

if  the  Jacobian  changes  its  sign,  and  by  n,  if  both  momenta  change  their  sign  during  the  same  time 
step. 


D.  Average  of  the  final  amplitudes  over  the  initial  conditions 
The  exact  quantum  mechanical  wave  function  has  the  asymptotics: 
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where  the  S-matrix  elements  S  %  y  ^ ,  for  vj  *  v’j'  are  given  by  the  relation: 
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By  analogy  with  Eq.  (154)  we  evaluate  the  semiclassical  S-matrix  elements  by  means  of  the  relation: 
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To  avoid  of  the  singularities  in  the  turning  points  we  first  represent  this  matrix  element  as  the  two- 
dimensional  integral  using  the  relation: 
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(The  relation  turns  into  the  identity  if  we  deal  with  the  exact  matrix  element  which  is  independent  of 
R.)  Substituting  Eq.  (155)  into  Eq.  (156)  we  find 
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with 


(^(t,r0,t>)  =  P1(t0,r0,\))Qj(t,rg\)),  Q^t.r^u)  s  P1(t0,r(yo)Q11(t,r{))-u)  .  (160) 


and  the  notations  [t.'u]  and  [t,x,D]  stand  for  ( rQ(T),u)  and  (t,  ^(t),;)).  We  directly  propagate 
-0  _1 

Q  j  and  along  the  trajectory  to  avoid  possible  complications  coming  from  the  singularities  in 
Eq.  (143)  which  appear  if  the  trajectory  starts  from  the  turning  point. 

3.  P.9t?mial  Energy  gyring 

Two  candidates  of  collisional  excitation  of  vibrational  levels  were  studied,  (1)  in 
O(3P)+HF(v=0j=0)  collisions  and,  (2)  in  Na(3s2S)+H2(v=0j=0)  collisions,  both  of  which  have 

previously  existing  potential  energy  surface  information.  The  potential  energy  surface  for  the  OHF 
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system  is  based  on  an  earlier  UHF  calculation2  consisting  of  a  sum  over  pairs  plus  three-body 
correcting  analytical  fit  (surface  1  of  Ref.  2)  whereas  the  NaH2  collisions  were  studied  on  a  DIM 

potential  energy  surface.6  A  LEPS  fit  to  MCSCF  results  for  OHF  were  also  generated  in  Ref.2  to 
enable  the  dynamical  study  of  the  H  abstraction  reaction,  but  this  needs  to  be  employed  only  at 
higher  energies  than  presently  reported.  We  refrain  from  including  the  analytical  and  numerical 
details  regarding  the  potential  surfaces  for  brevity  since  they  are  available  elsewhere. 

A  central  feature  of  the  OHF  potential  surface  that  is  relevant  for  vibrational  excitation 
dynamics  is  the  expected  concentration  of  most  of  the  excitation  (at  lower  energies)  to  a  narrow 
range  of  IOS-angles,  close  to  the  approach  of  O  atom  from  the  H-atom  end  of  the  HF  molecule;  in 
this  direction,  the  potential  surface  allows  a  reactive  encounter  to  occur  at  higher  energies.  This  • 
will  be  borne  out  by  the  dynamics  results  below.  In  the  case  of  NaH2  collisions,  only  the 

electronically  adiabatic  ground  surface  (obtained  from  diagonalizing  the  eight  diabatic  DIM  surfaces 
of  Ref.  6)  is  employed  in  the  electronically  elastic  results  to  be  presented  below  for  comparing 
QIOS  and  RSVE  results. 

4.  Calculational  Details 

The  RSVE  dynamical  calculations  using  Eqs.  (129)  are  straightforward  once  the  potential 
matrix  in  the  vibronic  basis  is  set  up.  This  is  done  by  diagonalizing  the  asymptotic  vibrational 
Schroedinger  equation  in  a  harmonic  basis  set  or  employing  one  of  many  standard  numerical 
algorithms1  to  obtain  the  target  states.  The  atom-diatom  potential  matrix  is  readily  computed  along 
the  trajectory  since  it  is  expanded  in  these  target  states  which  have  amplitudes  that  vary  with  time 
according  to  Eq.  (129),  thereby  leading  to  a  time-dependent  potential  obeying  Eq.  (128).  The 
amplitude  initial  conditions  are  chosen  to  correspond  to  unit  probability  in  the  diatomic  v=0  state. 
The  separation  R  is  chosen  large  and  the  asymptotic  negative  relative  momentum  is  defined  by  the 
channel  potential  in  the  v=0  channel  and  collision  energy.  The  final  amplitudes  in  other  state 
channels  resulting  from  solving  the  dynamical  equations  yield  the  inelastic  transition  probabilities. 
The  computation  of  cross  sections  involves  integration  over  many  trajectories  that  represent 
varying  IOS-angles  and  impact  parameters  for  each  given  initial  rovibrational  state.  All  the 
preliminary  computations  reported  here  are  for  initial  v=j=0  for  the  diatomic  and  the  study 
examines  dependence  on  collision  energy  and  total  angular  momentum  (semiclassically  defined  by 
the  impact  parameter). 

The  SIOS  calculations  start  with  the  generation  of  a  bundle  of  classical  trajectories 
determined  by  the  collision  energy  and  total  angular  momentum.  Only  20  trajectories  are  employed 
in  the  present  results.  The  main  task  involves  building  the  semiclassical  eikonal  wavefunction 
along  these  trajectories  giving  special  attention  to  the  semiclassical  phase  changes  at  caustics  (such 
as  classical  turning  points).  The  latter  are  determined  by  the  semiclassical  Bom  Oppenheimer 
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approximation  (see  previous  section  and  Appendix  2)  applied  to  the  present  problem.  The  post¬ 
collision  eikonal  wavefiinction  is  then  projected  onto  selected  final  states  to  directly  determine  the 
S-matrix  elements  for  inelastic  transitions.  The  projection  process  employs  interpolation  of  the 
asymptotic  eikonal  wavefunction  to  result  in  a  finely  grided  numerical  quadrature  for  S-matrix 
amplitudes  since  the  inelastic  trajectories  end  in  a  nonuniform  .grid  even  when  the  initial  conditions 
start  from  a  uniform  grid.  Cross  section  computations  involve  integration  of  these  amplitudes  over 
IOS-angles  and  impact  parameters  and  introduce  further  quantal  interference  effects. 

5.  Results  and  Discussion 

Figure  6  shows  that  excellent  agreement  with  an  exact  limit  QIOS  benchmark  is  obtainable 
from  the  RSVE  computations.  This  success  of  RSVE  in  gas  phase  is  consistent  with  the  similar 
success  that  the  coupled  states  ansatz  has.  had  in  treating  rotationally-electronically  inelastic 
gas/surface  encounters  involving  up  to  282  channels.7  The  most  exciting  aspect  of  Fig.  6  is 
however  the  promise  shown  there  by  the  SIOS  method,  which  is  seen  to  capture  the  essential 
physics  of  this  collision;  the  results  shown  are  from  SIOS  calculations  that  employ  semiclassical 
wavefunctions  based  on  a  meager  20  trajectories  and  are  not  yet  numerically  the  best  they  can  be. 

O+HF  Vibrational  Excitation 
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Figure  6.  Comparison  of  angle-dependent  accurate  quantum-IOS  probabilities  with  those 
obtained  with  the  new  semiclassical  method  for  v=0  — » v=l  vibrational  excitation 

in  collisions  between  0(3p)  atoms  and  HF  molecules  at  a  translational  energy  of 
3  eV.  Both  total  and  internal  angular  momenta  are  zero. 
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Figures  7  and  8  show  dependence  on  total  angular  momentum  for  the  two  examples  of 
vibrational  excitation  studied  here.  Fig.  7  contains  all  the  data  available  at  the  time  of  this  report. 

(Detailed  calculations  are  in  progress  to  make  a  full  assessment  of  the  methodologies  in  different  L 
and  Etr  regimes.)  The  NaH2  computations  were  done  to  see  if  the  quality  of  agreement  is 

maintained  for  a  completely  different  molecular  system  and  Fig.  8  is  encouraging  (note  that  H2 
being  a  homonuclear  diatomic  molecule  is  fully  studied  by  sampling  half  the  lOS-angle  space 
compared  to  HF). 


O+HF  Vibrational  Excitation 


Figure  7.  Comparison  of  angle-dependence  of  accurate  quantum-IOS  opacities  with 
those  obtained  with  the  new  semiclassical  method  for  O(-^P)  atoms  colliding 
with  ground-state  HF  molecules.  Curves  correspond  to  different  values  of 
orbital  angular  momentum.  Translational  energy  is  3.0  eV. 
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Na+H2  Vibrational  Excitation 


Figure  8.  Comparison  of  angle-dependent  opacities  from  IOS  quantum  calculations  with 
those  obtained  with  the  new  semiclassical  method  for  (3s^S)  Na  atoms  colliding 
with  ground-state  H2  molecules.  Curves  correspond  to  different  values  of  orbital 
angular  momentum.  Total  energy  is  2.04  eV. 

6.  Concluding  Remarks 

We  have  produced  a  systematic  numerical  benchmark  for  collisional  excitation  of 
vibrational  levels  employing  semiclassical  methods.  Two  semiclassical  computational  routes,  (1) 
using  a  coupled  quantal  states  ansatz  to  describe  vibrational  states  quantum  mechanically  and 
translational  motion  classically,  within  the  SCEM  framework,  and  (2)  using  semiclassical 
wavefunctions  generated  from  classical  trajectory  results  using  a  newly  developed  (semiclassical 
Bom-Oppenheimer  approximation)  procedure  were  tested  against  exact  quantum  limit  results.  Both 
the  semiclassical  methods  were  successful  in  obtaining  reliable  numbers  for  state-to-state 
vibrational  excitation  probabilities.  The  use  of  method  (2),  termed  SIOS  herein,  is  indicated  as  an 
exciting  prospect:  such  a  technique  is  to  be  welcomed  as  a  computational  route  allowing  the 
transformation  of  information  from  classical  trajectories  to  state-to-state  transition  amplitudes  at  the 
semiclassical  level.  The  novel  feature  of  the  SIOS  scheme  employed  here  compared  to  rigorous  S- 
matrix  theory  of  Miller8  (which  was  difficult  to  apply  in  practical  problems  and  soon  went  out  of 
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use)  is  that  there  are  no  nonlinear  root  trajectory  searches  to  be  performed  here;  instead,  amplitudes 
are  obtained  by  projecting  the  eikonal  scattering  wavefunction  on  to  specific  final  states.  The 
promise  of  such  a  method  was  demonstrated  earlier  in  the  Franck-Condon  region  for 
photodissociation,7 8 9  but  the  present  validation  for  vibrational  excitation  processes  significantly 
extends  the  validity  regime,  by  demonstrating  that  useful  quality  persists  in  outgoing  eikonal 
wavefunctions  in  the  asymptotic  region  provided  valid  connection  formulae  are  employed. 
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The  Quenching  of  Na(32P)  by  H2:  A  Quantal  IOS  Calculation  of 
Electronic-to*Vibrational  Energy  Transfer 

i,  Introduction 

The  quenching  of  electronically  excited  Na  by  H2  is  a  prototype  for  elec tronic-to- vibrational 
(E-V)  energy  transfer  in  molecular  collisions.  E-V  energy  transfer  involving  alkali  atoms  has  been 
widely  studied  both  experimentally  and  theoretically.1'14  The  theoretical  interpretation  of  the  Na- 
H2  quenching  process  has  focused  upon  the  two  lowest  states  of  2A'  symmetry;6’10’13’15'17  for 
C2v  geometries  these  are  the  12Aj  and  12B2  states.  The  2B2  state  has  a  well;  it  dissociates  to 
Na(3p2P)  +  H2,  and  it  exhibits  a  conical  intersection  with  the  12Aj  surface  that  dissociates  to 
ground-electronic- state  reagents.  For  less  symmetric  geometries,  both  these  states  become  2A\ 
Previous  dynamical  studies  indicate  that  these  two  electronic  states  can  qualitatively  describe  the 
quenching  process.5,10,13  A  quantitative  description  of  this  process  requires  both  potential  energy 
surfaces  for  the  electronic  states  of  importance  and  the  nonadiabatic  coupling  terms  between  these 
electronic  states.  Ab  initio  restricted  Hartree-Fock  calculations  have  been  reported  for  the  potential 
energies  for  the  two  lowest 2 A'  states  for  about  100  nuclear  arrangements;  however,  no  ab  initio 
calculations  of  the  nonadiabatic  coupling  terms  have  been  reported,  and  the  potential  energy 
calculations  including  electron  correlation  (by  the  coupled-electron-pairs  approximation  -  CEPA) 
have  been  performed  only  for  geometries.6,10 

Recently  the  diatomics-in-molecules  (DIM)  method18,19  has  been  applied  to  calculate  the 
three  lowest-energy  2A’  potential  energy  surfaces  of  NaHj.13,17  The  two  lowest-energy  surfaces 

agree  well  with  the  available  ab  initio  calculations  of  Botschwina  and  Meyer.6,10  A  major 
advantage  of  using  the  DIM  formalism  is  that  it  provides  a  global  representation  of  the  potential 
energy  surfaces  and  the  couplings  between  them. 

There  has  been  a  large  amount  of  work  on  developing  methods  for  the  quantum  mechanical 
treatment  of  electronic  transitions  in  atom-molecule  collisions.20'32  Zimmermann  and  George21 
compared  the  use  of  diabatic  and  adiabatic  representations  for  atom-vibrator  collisions.  In  the 
adiabatic  representation  coupling  between  different  electronic  states  introduces  first-derivative 
terms  in  the  coupled-channel  equations.  Zimmermann  and  George  solved  the  adiabatic  coupled- 
channel  equations  by  transforming  the  N  second-order  differential  equations  into  a  set  of  2N  first- 
order  differential  equations.  Another  alternative,  employed  by  Baer  and  coworkers,23'25  is  to 
obtain  coupled-channel  equations  that  could  be  solved  using  the  efficient  computational  algorithms 
developed  for  second-order  differential  equations  without  first-derivative  terms.  In  Baer's 
approach  a  transformation  is  made  from  the  adiabatic  representation  to  a  diabatic  representation  in 
which  all  coupling  arises  through  off-diagonal  elements  of  the  potential  matrix.  In  general  such  a 
transformation  does  not  exist,33  but  it  does  exist  if  the  nonadiabatic  couplings  are  approximated  by 
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the  Preston-Tully19*34  approximation,  as  used  in  the  DIM  method.  Baer's  approach  has  been 
applied  to  collinear  atom-diatomic  collisions,  and  more  recently  Rebentrost  and  Lester26*27  have 
used  a  similar  diabatic  transformation  to  treat  the  nonreactive  F  +  H2  problem  in  three  dimensions. 

The  previous  quantum  mechanical  calculation  on  Na  +  H2  by  McGuire  and  Bellum5  used 
the  two  lowest-energy  adiabatic  potential  energy  surfaces  as  an  approximation  to  the  diagonal 
diabatic  potential  energy  surfaces.  The  off-diagonal  terms  of  the  diabatic  potential  were  modelled 
as  gaussian  terms  centered  on  the  location  of  the  minimum  of  the  energy  difference  of  the  adiabatic 
surfaces.  Cross  sections  were  calculated  within  the  rotational  infinite-order-sudden  (IOS) 
approximation;35*36  however,  the  integral  over  the  IOS  angle  was  approximated  using  only  the 
perpendicular- bisector  approach  of  the  Na  to  the  H2  molecule. 

In  the  present  calculation  we  describe  a  new  method  for  calculating  electronic  transition 
cross  sections  using  adiabatic  potential  energy  surfaces  and  the  nonadiabatic  coupling  terms  as 
input.  In  this  approach  we  use  a  mixed  adiabatic-diabatic  representation  in  which  the  motion  in  the 
diatomic  intemuclear  distance  is  treated  diabatically  and  the  motion  in  the  Na  to  center-of-mass  of 
H2  distance  is  treated  adiabatically.  We  use  the  rotational  IOS  approximation  to  decouple  the 

scattering  dynamics  in  the  orientation  angle.  The  R  matrix  propagation  method30"32,37*42  is  used 
to  solve  the  coupled  equations  describing  the  scattering  process  in  the  mixed  representation  for 
each  value  of  the  orientation  angle,  and  the  calculations  are  converged  with  respect  to  the 
numerical  integral  over  the  orientation  angle. 

In  Subsection  2  we  present  the  details  of  the  new  theoretical  method,  and  in  Subsection  3 
we  present  details  of  the  computational  procedure.  Subsection  4  contains  results  of  the 
calculations,  and  Subssection  5  presents  a  discussion. 


2,  Theory 

We  consider  the  collision  of  an  atom  A  with  a  diatomic  molecule  BC,  where  R  is  the 
distance  from  A  to  the  center  of  mass  of  BC,  r  is  the  BC  intemuclear  distance,  and  y  is  the  angle 
between  the  R  and  r  vectors.  Within  the  rotational  IOS  approximation  the  total  Hamiltonian  is 
given  by36 


R  + 


h2 1  (I  +  1)  . 

2  + 

2pR 


(161) 


where  |i  is  the  reduced  mass  for  the  R  motion,  I  and  j  are  particular  values  of  the  orbital  and 

rotational  quantum  numbers,  respectively,  and  the  Hamiltonian  for  vibronic  coordinates  (x  and  r 
where  x  is  the  collection  of  electronic  coordinates)  H^  is  given  by 
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+  H  , 
el 


(162) 


Here  jxBC  is  the  reduced  mass  for  the  r  motion  and  Hel  is  the  electronic  Hamiltonian.  The 
approximate  Bom-Oppenheimer  electronically  adiabatic  eigenfunctions  <{>n  for  electronic  state  n 
satisfy  the  matrix  eigenvalue  equations 
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where  Vn  (r,R,y)  is  the  electronically  adiabatic  potential  energy  surface  for  electronic  state  n,  (t) 
denotes  dependence  on  all  coordinates  (x,r,R,y),’and  we  use  a  bra-ket  notation  to  denote  integrals  - 
-  the  subscript  denotes  the  variables  integrated  out.  Within  the  rotational  IOS  formalism  the  orbital 
angular  momenta  is  conserved  and  the  rotational  angular  momentum  is  not  coupled  with  the 
electronic  angular  momenta. 

A  mixed  basis  for  the  electronic  degrees  of  freedom  is  obtained  as  follows.  For  fixed  R 
and  y  a  transformation  is  made  to  a  representation  {<j>^ }  which  is  P-diabatic31*43  with  respect  to 

the  r  motion  but  adiabatic  with  respect  to  R: 
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where  n,^  Note  that  unless  nj^  is  infinite,  0n>  is  actually  a  function  of  r  as  well  as  x, 

R,  and  y.  This  dependence  is,  however,  effectively  removed  in  the 
nmax  'fold  subspace  of  retained  adiabatic  functions  by  the  condition  of  P-diabaticity,  namely. 


n,n'=l 


’•",nmax 


(166) 


The  transformation  from  the  adiabatic  to  the  mixed  basis  is  given  by 


3  max  * 

aF  un'n(r’R’Y)  =  "  1  f^n"(r-^Yl  u„-n(r’R’Y) 
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where 


C„(r’R'?>  -  (<M  £!♦*«>), 


(168) 


We  require  u  to  be  an  orthogonal  matrix  so  that 
<  *^(T)  |  C (t)  >x=8n.n 


(169) 


Note  that  Eq.  (167)  is  analogous  to  a  well  known  equation  used  to  define  P-diabatic  bases  in  atom- 
atom  collisions.3 1‘33,44 

The  internal  wavefunctions  ^(t)  are  expanded  in  the  mixed  basis 


n _ _  M_ 

max  n 


Vm(x)=r1  I  I  <|>f(T)pJn  (r)CJ„  (R,y) 
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m=l,...,M 


(170) 


where  the  total  number  of  channel  is  given  by 


M  =  I  (M„  +  l) 
n=l 

and  pj,v(r)  is  a  vibrational  basis  function  for  electronic  state  n  satisfying 

(  r-1  p-*  ,  ,(r)  r-1  p„v(r) )  =  S  .  , 

v  r  n  v  '  rnv'  'if.  nvnv 


(171) 


(172) 


where  SnVnv=5v.y  for  n'=n.  Note  that  integration  over  r  includes  an  r2  volume  factor.  The 
internal  eigenvectors  are  defined  by 


m,m'=l,...,M 


(173) 


<<■«  >,,r  =  5 


m,m'=l,...,M 


(174) 


which  yields  the  eigenvalue  equation 


hVy)  CJ(R,y)  =  CJ(R,Y)  EJ(R,y) 
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(175) 


where  Ei  is  a  diagonal  matrix,  the  matrix  elements  of  the  internal  Hamiltonian  are  given  by 

HJ„.„  „„  <R.Y>=  T*  Pj„v  (r)| 
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and  the  potential  matrix  V^.  is  given  by 


V'V.R.y)  =  uT(r,R,y)  V*\r,R,y)  u(r,R,y) 


(176) 


(177) 


The  total  scattering  wavefunction  for  an  initial  state  cOq,  a  fixed  value  of  the  conserved 
quantity  I,  and  a  typical  or  average  value  j  of  the  rotational  quantum  number  is  expanded  in 
eigenvectors  of  the  internal  Hamiltonian 


fil.«)  =  n'1£vi.w  xL  (R.Y) 
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The  rotational-IOS  close-coupling  equations  for  a  fixed  y  arc  given  by 

fdJdrrVm'  «[Hlw-E]>r,io«-0 

Substituting  Eqs  (161)  and  (178)  into  Eq.  (179)  gives 
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The  coupling  terms  arc  defined  by 
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These  expression  can  be  evaluated  by  using  Eq.  (170);  however,  this  requires  calculations  of 
derivatives  of  (^(R.y)  and  matrix  elements  of  the  derivative  operators  in  the  0^  basis. 

The  radial  wavefunctions  are  subject  to  the  standard  boundary  conditions36 


J1 


-mm 


(R.Y) 


R— >0 


->  0 


(183) 


and 


J1 


rmm„(R.Y)  R-_k-^5mmoexp[-i(kjmR  -  ln/2); 


-SJL  W  cxp[i(kjmR  -  !  n/2)] 


(184) 


where  the  asymptotic  wavenumber  is  defined  by 


(185) 


Opacity  functions  can  be  defined  for  each  value  of  j  and  I 


^  ,  =  —  j  dy  sin  y 

m0^m  2jq1 


S.  -SJ‘  (y) 

1  m  m  q  m  m  q 


(186) 


and  the  total  cross  section  for  transition  from  initial  state  j0,  to  final  level  m'  summed  over  all 
final  rotational  quantum  numbers  j'  is 


crj_  .  =  -t—  I  (21  +  l)  p*.1 


m0->m  k2 

jm, 


m  0~>m 


(187) 


Notice  that  this  result  is  independent  of  j0  in  the  rotational  IOS,  but  it  does  depend  slightly  on  j, 
which  is  a  parameter  of  the  approximation  scheme. 
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3.  Calculational  details 


A.  Vibrational  basis  functions 

For  the  Na  +  H2  system,  the  electronic  states  considered  here  should  correlate  with  H2  in 
its  ground  electronic  state  when  the  Na  to  H2  distance  is  large.  However,  within  the  DIM  model 
used  here  the  second  electronic  state  corresponds  to  Na  (3p2P)  +  H^Zg)  for  r  less  than  about  3.1 
Oq,  but  for  r  greater  than  this  value  it  crosses  with  Na(3s2S)  +  H2(3IU).  At  the  crossing  point  the 
energy  is  about  5eV  above  the  zero  of  energy  at  Na(3s2S)  +  H^Zg)  with  r  =  req,  but  the  potential 
then  drops  to  about  4.74  eV  as  r  increases.  Therefore,  the  vibrational  states  are  computed 
individually  for  the  two  electronic  states.  Within  the  DIM  formalism,  all  electronic  coupling 
vanishes  for  Na  infinitely  separated  from  H2  and  the  asymptotic  diatomic  potential  are  simply  the 

adaibatic  potentials  for  the  vibrating,  rotating  diatomic  molecule.  The  vibrational  basis  function  are 
defined  by 


ft2j(j  +  1) 

+  -;  ■  r  +  VHHn(r) 

2d  QCr2  HH’n 


Pnv(r)  ^nv  Pnv(r) 


(188) 


where  VHH  n(r)  is  the  asymptotic  H2  vibrational  potential  for  the  electronic  state  n.  The  energy 
levels  £[,  are  obtained  numerically  by  the  Cooley-Cashion  method.45  All  integrals  over  r  are  done 
by  Gauss-Hermite  quadrature  using  50-100  points.  The  quadrature  points  are  defined  by 
ri~req  +  x/  a  where  Xj  are  the  Gauss-Hermite  quadrature  nodes,  and  a  is  the  range  parameter 
for  the  ground  state  H2  vibrational  potential  taken  to  be  4.57082  a^.  The  numerical  values  of  the 
vibrational  basis  functions  are  stored  on  the  grid  of  quadrature  points  { rk }  j  to  be  used  in 
subsequent  numerical  integrals,  e.g.  Eq.  (176). 

B.  P-diabatic  transformation 

Equation  (167)  for  the  transformation  matrix  u  is  solved  by  the  Magnus  method.32  The 
transformation  matrix  is  only  needed  at  the  Gaussian  quadrature  points  and  the  Magnus 
approximation  yields 

U(rk+1,R,Y)  =  exp  {- (rk+l  -  rk)  +rk),R.Y]}»(rk.R,Y)  g9 


The  choice  of  u  at  the  first  grid  point  is  arbitrary,  and  if  no  further  approximations  are  made  the 
final  results  of  the  calculations  are  independent  of  this  choice.  For  convenience  we  set 
unn.(r0,R,Y)=6nn-  for  all  R,y  and  find  u  for  r<r0  and  r>r0  by  inward  and  outward  application  of  Eq. 
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(185),  respectively.  r0  is  chosen  to  be  a  grid  point  near  re.  The  exponential  of  a  matrix  is 
evaluated  by  the  method  of  ref.  46. 


C.  R-matrix  propagation 

Equation  (180)  is  not  solved  direcdy  for  the  radial  wavefunctions;  instead  we  use  the  R- 
matrix  propagation  method30'32,37"42  to  obtain  a  global  R  matrix  for  relating  the  ratio  of  the  radial 
wavefunction  and  its  derivative  at  large  R  values  to  the  ratio  at  small  R  values.  The  S  matrix  is 
then  obtained  from  the  R  matrix.  The  method  for  treating  electronic  transitions  in  systems  with  one 
intemuclear  degree  of  freedom  by  the  R-matrix  propagation  method  when  the  input  is  in  an 
adiabatic  representation  has  been  presented  previously.31,32  This  method  is  extended  here  to  treat 
the  case  of  several  nuclear  degrees  of  freedom  where  the  representation  is  adiabatic  in  the 
propagation  coordinate  and  diabatic  in  other  coordinates,  as  in  Subsection  2.  In  the  R-matrix 
propagation  method  we  use  internal  basis  functions  \|^(x,r,RVy)  that  are  independent  of  the  radial 
coordinate  R  within  each  sector  (i)  centered  at  Rl.  Propagation  of  the  global  R  matrix  across  a 
sector  is  accomplished  in  terms  of  the  eigenvalues  E^fR'.y)  of  the  internal  Hamiltonian.  The 
coupling  between  the  internal  states  arises  at  the  sector  boundaries  in  transforming  from  the  internal 
basis  functions  in  one  sector  to  those  in  the  adjacent  sector.  It  is  the  definition  of  this  sector 
transformation  matrix  which  is  the  essential  aspect  of  the  scheme. 

We  will  consider  two  ways  to  calculate  the  sector  transformation  matrix  with  no  further 
approximations,  i.e.,  so  that  the  results  are  independent  of  the  initial  condition  used  in  integrating 
Eq.  (189)  and  are  equivalent  to  what  would  be  obtained  by  numerical  integration  of  the  Eqs.  (180). 
Then  we  will  consider  a  fixed-diabatic-states  approximation  in  which  the  P-diabaticity  in  the  r 
coordinate  of  the  basis  functions  of  Subsection  2  is  interpreted  as  if 


JL 

dr 


(x)  =0 


(190) 


Assuming  Eq.  (190)  is  mathematically  equivalent  to  assuming  that  the  electronic  states  included  in 
the  expansion  are  a  complete  set  since  Eq.  (166)  shows  that  the  derivative  of  Eq.  (190)  is 
orthogonal  to  the  other  electronic  functions  included  in  the  electronic  basis.  The  fixed-diabatic- 
states  approximation  greatly  simplifies  the  numerical  work;  we  shall  test  its  accuracy  by  comparing 
the  final  cross  sections  to  those  computed  without  this  approximation. 

First  consider  the  accurate  calculation  of  the  sector  transformation  matrix.  The 
transformation  matrix  from  sector  i  to  sector  i+1  is  defined  by 


T!n’m1(Y)  =  (  ^mXx,x,R\y)  |  ^m(x,r,R1+\y)  )x  , 


(191) 
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This  can  be  approximated  by31,32 
i,i+l 


T1’!  (Y)  =  exp  {(r‘+1 -Rl)  FJ  ,  T4(r‘  +  R1+1),yl} 

m  m  "  K  L  m  m.2  ’*  JJ 


(192) 


where  F^m.m(R,Y)  is  defined  in  Eq.  (181).  This  expression  requires  the  evaluation  of  numerical 
first  derivatives  and  to  avoid  having  to  compute  the  numerical  derivatives  of  u(r,R,y)  and  C(R,y) 
we  use  an  alternative  method  to  obtain  the  sector  transformation  matrix.  Substituting  Eq.  (170) 
into  Eq.  (192)  yields 


i,i  +  l  "max  ^n'  "max  Mn  ;  •  i,i+l 

"*"m'm^  =  ?  ^  ^  ^  Cn'v  „m^R  ^n'v  „nv  ^  C*  (R1  y)  (193) 

n  =1  v ^,=0  n=l  v  n=0  n  n'  n  nv  , m,'‘  ’ r 


where  the  overlap  matrix  is  defined  by 

i  ,i  +1  00  i  dA  i  dA  i  +1 

°n'v  nv  (Y)  =  ldrPv  (r)(<t>n,(x,r,R  ,y)  <J>n(x,r,R  ,y)  >x  pJy  (r) 
"  0  " 


n'  n 


(194) 


Using  Eq.  (165)  gives  our  final  expression  for  the  overlap  matrix 


i  ,i+l 


»  i*  '  *  *  -  I  |  |  "T'X 

°„-v  nv  (y)  =  ldrpJv  (r)  I  un„n,(r,R,y)  t;„n„,(r,y) 
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n  v  ,nv_ 
n  n 


n  n"n'"  n  n 


i,i+l 
n  ”n' 


%-"n<r’R,+‘r)P,»„<r) 


where 


t'n.'n+l(r.r)  »  <  <t>^'(x,r,R1,Y)  |  o“(x.r,R‘+1,y)  >„ 

=  «p{(Ri+1-R‘)  F^n[r,^(R'  +  R'+1).y]} 

-  (*>>  I  ♦?<*>)„ 


and 


(195) 

(196) 

(197) 


(198) 


It  is  important  to  note  that  the  overlap  matrix  in  Eq.  (195)  is  not  necessarily  orthogonal. 
Orthogonal  overlap  matrices  arc  guaranteed  only  when  the  primitve  basis  functions  in  this  case) 
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do  not  change  from  sector  to  sector.  However,  because  of  the  diabatic  transformation  in  each 
sector  the  overlap  matrix  is  orthogonal  only  under  special  conditions.  Consider  the  matrix  product 


n  N  „  .  .  , 

max  n  i,i+l 


?.  ^  n°n”v  „n'v  °n"v  , 

n  =1  v  =0  n  n  n 
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n  M  .. 

max  oo  oo  .  .  n 


=  I  /  dr  /  dr  ’  pJv  (r)pJv  (r’)  I  pJv  (r)pJv  (r1) 
n”=l  0  0  n‘ 


n  v  =0  n" 
n 


♦?•<*>  >*<♦*«>  I  ♦*<i+1>>i  (199) 


dA  dA 

where  <()n  (i)  denotes  (f>n  (x,r,R\y)*  The  overlap  matrix  is  orthogonal  only  if  the  matrix  defined  in 
Eq.  (199)  is  the  unit  matrix.  This  can  happen  only  if  the  vibrational  states  within  each  electronic 
manifold  are  a  complete  set,  i.e.. 


P  v  (r)p  y  (r’)  =  S(r  -  r’) 

n”  n" 


(200) 


and  the  electronic  states  are  also  a  complete  set,  i.e., 
"r  I  .  dA...  .  ,  ,  dA,.,  .  , 

I  ♦n  W><  <t>n  (»)  I  =1 

n=l 


(201) 


For  the  calculations  performed  here  these  conditions  are  not  generally  valid,  and  we  therefore  use 
the  correct  equations  for  propagation  of  the  sector  R  matrices.31,47 

If  the  fixed-diabatic- states  approximation  is  employed,  Eq.  (190),  the  matrix  elements 

i,i+l 

1  n  n  ^r,Y^  become  independent  of  r  and  can  be  evaluated  at  one  arbitrary  value  of  r.  For  the 
chose  of  r*=r0,  the  expression  for  the  sector  transformation  matrix  simplifies  to 


i,i+l  M  ;  .  n  max  • 

Ttn'm(Y1  *  ?  Cp'm'(R  ’Y)  5v(p’)v(p)  J  U n'n(p')(rcr R  ’Y) 
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(202) 
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4.  Results 

21.A. 

In  Fig.  9,  contours  of  the  energy  for  the  two  lowest  adiabatic  potentials  Vn  (r,R,y)  are 

shown  as  a  function  of  the  H-H  intemuclear  distance  r  and  the  distance  R  from  the  Na  atom  to  the 

center  of  mass  of  H2,  at  three  values  of  the  angle  y.  In  Fig.  10,  the  two  lowest  adiabatic  potential 

energies  along  cuts  in  r  are  show  at  three  different  values  of  R  and  for  the  same  three  values  of  y. 

Evidence  of  an  avoided  crossing  between  the  second  state  and  third  state  are  evident  for  the 

potentials  at  R=5  aQ.  The  nonadiabatic  coupling  element  in  r  between  these  two  lowest  states 

f^(r,R,y)  is  shown  in  Fig.  11  for  the  three  values  of  R  at  the  three  values  of  y.  The  coupling 

exhibits  one  major  peak  which  shifts  to  higher  values  of  r  for  increasing  values  of  R.  The  mixed 

dA 

adiabatic-diabatic  potential  matrix  Vnn  (r,R,y)  obtaining  from  the  P-diabatic  transformation  is 
displayed  in  Fig.  12.  The  diagonal  elements  are  close  to  the  adiabatic  curves  shown  in  Fig.  20 
especially  near  the  minima  in  the  potentials.  The  two  diagonal  curves  tend  to  cross  near  values  of  r 
which  displayed  maximum  in  the  nonadiabatic  coupling.  Figure  13  shows  the  eigenvalues  of  the 
internal  Hamiltonian  EJm(R,y)  as  a  function  of  the  Na  to  H2  distance  for  a  fixed  value  of  y.  These 

curves  exhibit  multiple  avoided  crossing  indicative  of  a  system  with  strong  coupling.  Figure  14 
shows  the  first  preliminary  results  for  scattering  of  Na(3p2P)  from  in  its  ground  vibrational 

state  at  a  total  energy  of  3  eV.  Probabilities  are  shown  for  only  the  4  final  states  with  the  largest 
probabilities.  The  elastic  channel  is  dominant,  but  the  v=2  state  is  shown  to  be  much  larger  than  all 
other  states  except  for  v=0  and  4  at  angles  near  30°.  This  is  qualitatively  in  agreement  with  the 
experimental  results  of  Hertel  et.  al.8,9 
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Figure  9.  Potential  energy  contours  for  the  two  lowest  adiabatic  potential  energy  surface  (r,R,y)  of  NaH2  as  functions  of  the  H2  distance  r  and 

Na-to-H2  distance  R,  for  fixed  angle  y  between  the  r  and  R  vectors.  The  left,  center,  and  right  columns  of  plots  are  for  y  =  1 5, 45,  and  75 
degrees,  respectively  (90  degrees  corresponds  to  perpendicular  approach  of  Na  to  H2).  The  top  row  displays  the  lowest  adiabatic 

potential.  The  contours  are  evenly  spaced  at  0  5  eV;  the  zero  of  energy  is  Na(3  S)  infinitely  separated  from  H2  at  its  equilibrium  geometry. 


Q  A 

Figure  10.  Adiabatic  potential  energy  curves  Vn  (r,R,y)  for  the  two  lowest  energy  states  of 
NaH2.  The  curves  are  plotted  as  a  function  of  the  H2  distance  r  for  fixed  Na-to-H2 

distance  R  and  for  fixed  angle  y  between  the  r  and  R  vectors.  The  top,  middle,  and 
bottom  rows  of  plots  are  for  R  =  3,  4,  and  5  a^,  respectively.  The  left,  center,  and 

right  columns  of  plots  are  for  y  =  15,  45,  and  75  degrees,  respectively  (90  degrees 
corresponds  to  perpendicular  approach  of  Na  to  H2). 
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Figure  1 1.  Same  as  figure  10,  expect  for  the  nonadiabatic  coupling  terms  f^(r,R,y)  between  the 
two  lowest  energy  states  of  NaH2. 
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dA 

Figure  12.  Same  as  figure  12,  expect  for  the  mixed  adiabatic-diabatic  potential  curves  Vnn.(r,R,y) 
obtained  from  the  two  lowest  adiabatic  states  of  NaH2  ^  their  nonadiabatic  coupling 
elements.  The  transformation  matrix  from  the  adiabatic  to  diabatic  representation  is 
arbitrarily  chosen  to  be  the  unit  matrix  at  r=1.4  a^. 
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Figure  13.  Eigenvalues  EJm(R,Y)  of  the  internal  Hamiltonian  in  the  mixed  adiabaic-diabatic 
electronic  representation  for  y  =  2.7°. 
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Figure  14.  Transition  probabilities as  a  function  of  the  IOS  angle  y  for  collisions  of 

Na(3p2P)  with  Hj  in  its  ground  vibrational  state  at  a  total  energy  of  3  eV.  The  solid 

curve  with  the  largest  magnitude  is  for  the  elastic  channel  --  production  of  Na(3p2P) 
and  H2  in  its  ground  vibrational  state.  The  long  dashed,  the  other  solid,  the  short 

dashed  curve  are  for  quenching  of  Na  to  the  3s2S  state  to  form  H2  v=0,  2,  and  4, 
respectively. 
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Condensed  Phase  Modelling  and  Computer  Experiments 

Introduction 

Condensed  phase  modelling  and  simulations  are  essential  to  HEDM  research.  The 
condensed  phase  will  be  the  relevant  phase  for  storage  and  initial  processing  of  HEDM's.  Our 
focus  is  on  treating  those  chemical  processes  involving  energetic  species  that  are  local  centers  of 
energy  storage  in  the  medium  (in  contrast  to  an  extended  storage  in  an  overall  metastable  solid 
matrix).  The  key  role  of  condensed  phase  modelling  and  simulations  in  our  research  can  be 
appreciated  if  the  system  is  imagined  to  be  divided  into  a  few  atoms  constituting  the  energetic 
species  and  the  rest  of  the  bulk  phase  material.  The  few  atoms  constitute  a  primary  region,  while 
the  bulk  phase  material  constitutes  a  heat  bath. 

The  present  three  year  program  employed  a  new  capability  involving  semiclassical 
methods.  The  new  method  treated  the  fully  correlated  dynamics  of  the  chemistry  in  the  primary 
zone  region.  If  this  region  is  only  subject  to  electronically  adiabatic  dynamics,  and  did  not  involve 
energy  transfer  between  quantized  nuclear  (rovibrational,  phonon)  states,  it  could  be  readily  treated 
by  electronically  adiabatic  classical  dynamics,  that  was  well-established  prior  to  our  research 
program.  When  the  key  to  energy  leakage  or  transfer  away  from  the  primary  zone  lies  in  quantum 
mechanical  effects,  involving  multidimensional  quantized  vibrations  or  possible  electronic 
inelasticity,  it  is  known  that  a. simple  classical  description  is  not  sufficient.  The  present  new 
semiclassical  methods  could  be  employed  to  treat  these  many-atom  energy  transfer  and  dynamics 
problems  since  quantum  mechanical  dynamical  techniques  are  not  computationally  feasible. 

One  role  of  our  simulation  work  is  to  define  parameters  to  describe  the  heatbath  dynamics 
properly  and  efficiently  via  few-body  heatbath  models;  the  latter  improve  the  efficiency  of 
performing  repeated  dynamics  calculations,  once  a  heatbath  is  characterized.  Another  important 
role  of  simulations  is  to  reveal  the  truth  about  what  happens  in  a  given  model;  since  simulations  are 
"computer  experiments",  they  often  reveal  unanticipated  phenomena.  Our  computer  experiments 
have  been  restricted  to  the  proposed  classical  Monte  Carlo  (MC)  and  Molecular  Dynamics  (MD) 
simulation  methods  from  the  outset  since  quantal  bath  effects  constitute  an  extensive  problem  by 
themselves.  Thus  the  studies  presented  in  this  report  employ  the  systems  helium  and  hydrogen  at 
high  pressures,  so  that  they  are  well  described  by  classical  simulations. 

The  following  subsection  presents  results  of  our  simulations  on  storing  excited  metastable 
helium  atoms  in  high  pressure  (GPa)  bulk  helium  liquid.  This  is  followed  by  a  subsection  on 
heatbath  modelling  based  on  the  methods  employed  in  this  research. 
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Monte  Carlo  Simulations  of  Helium  Bubble  States 

1.  Introduction 

It  is  known  from  studies  on  excited  helium  species1  produced  by  keV  electron 
bombardment  of  liquid  He  that  there  are  excited  atomic  and  molecular  species  trapped  in  physical 
bubbles,  the  latter  essentially  involving  a  vacuum  embedded  in  the  liquid  bulk.  The  nature  of  these 
bubble  states  have  been  discussed2  and  the  species  that  exist  within  them  have  been  identified  by 
spectroscopic  methods  and  subsequently  subjected  to  simple  model-based  interpretations.  The 
bubbles  are  stabilized  by  the  repulsive  interaction  of  the  Rydberg-like  excited  electron  with  bath  He 
atoms.  They  have  radii  in  the  range  5-20  A  depending  on  the  excited  species  and  the 
thermodynamic  state.  The  mechanism  that  supports  the  excited  species  in  a  bubble  is  closely  related 
to  that  for  the  case  of  electron-containing  bubbles  in  liquid  Helium.  The  repulsion  of  the  solvating 
He  atoms  of  the  liquid  by  the  excited  electron  of  the  atomic  or  molecular  species  is  considerably 
less  than  that  of  a  free  electron  since  the  nuclear  or  molecular  core  compensates  by  holding  on  to 
the  excited  electron.  Thus  the  radii  of  the  bubbles  in  the  electron  case  are  almost  twice  larger  than 
the  above.  A  related  phenomenon  of  recent  interest  has  been  that  of  noble  gas  bubbles  in  metals3 
where  there  have  been  a  few  atomistic  simulations.4  Atomistic  calculations  of  He*  bubbles  in 
liquid  He,  unlike  the  studies  of  noble-gas  bubbles  in  metals,  are  quite  limited.5 

The  previous  structural  analysis  of  the  bubble^  has  been  at  a  thermodynamic  state  point 
where  quantum  mechanical  effects  are  important.  The  calculations  were  carried  out  in  the  zero- 
temperature  limit,  from  first  principles,  using  a  variational  Jastrow  wave  function  and  the  integral 
equations  approach  for  fluid  mixtures.  It  may  be  asked  whether  these  quantal  liquid  effects  are 
always  necessary  for  the  bubble  formation.  When  the  thermodynamic  state  is  changed  by  raising 
the  temperature  or  pressure  or  varying  the  density,  regimes  will  appear  where  classical  dynamics  is 
appropriate  and  it  may  be  anticipated  from  simple  models2  of  the  classical  energetics  for  the  bubble 
formation  that  a  stable  bubble  will  be  formed  even  classically.  If  this  is  realistic,  it  means  that  a 
computer  simulation  methodology  based  on  classical  molecular  dynamics  and  Monte  Carlo 
techniques  will  provide  the  necessary  structural  and  dynamical  information  for  this  clasical  regime. 

The  classical  regime  occurs  at  room  temperature  in  a  very  high  pressure  range  (GPa 
regime).  In  this  regime  we  assess  the  influence  of  the  liquid  in  affecting  the  metastability  of  the 
electronically  excited  He-species.  Due  to  the  large  size  of  the  bubble  (of  diameter  5-10  A),  the 
trapped  excited  species  can  undergo  collisions  reminiscent  of  gas-surface  processes,  with  the  inner 
surface  of  the  bubble,  that  are  then  responsible  for  nonradiative  quenching.  Computer  simulations 
such  as  the  ones  explored  here  can  provide  the  data  (frequency  moments  of  the  heatbath  spectral 
density)  necessary  to  parameterize  a  Generalized  Langevin  Equation  (GLE)  model6  of  the 
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condensed  phase  dynamics  influencing  the  metastable's  quenching  process.  Our  plan  is  to  study 
the  electronically  inelastic  dynamics  of  the  trapped  metastable  species  by  the  Self-Consistent 
Eikonal  Method  (SCEM)7  in  a  merged  (SCEM+GLE)  stochastic  semiclassical  framework8. 

The  computer  simulations  also  help  explore  the  possible  phenomena  that  can  occur.  In  the 
present  system,  solvent  He  atoms  can  cluster  with  the  solute  He*  to  form  n-mers  depending  on 
pressure.  We  explore  the  n-values  of  the  Hen  clusters  occuring  in  the  bubble  states  under  the 

simulation  conditions  and  find  the  observed  average  <n>  as  the  pressure  is  increased.  At  very  low 
pressures  such  as  the  ones  corresponding  to  the  past  experiments  on  these  metastables,  primarily 
dimers  (n=2)  have  been  reported.  At  the  pressures  studied  here,  we  find  <n>  as  high  as  5-6  while 
the  bubble  state  character  of  the  (model)  system  is  still  in  tact.  The  <n>  value  is  also  dependent  on 
the  electronic  state  (a  or  c)  involved  and  thus  could  reveal  the  possible  importance  of  any  electronic 
nonadiabaticities  of  the  solute-solvent  interaction. 

Subsection  2  describes  the  interaction  model  employed  in  our  simulations,  and  Subsection 
3  gives  the  simulation  details  of  the  investigation.  The  concluding  Subsection  4  contains  the  main 
results  and  discussion. 

2.  Interaction  Potentials 

In  order  to  characterize  the  He*  bubble  in  liquid  He  for  various  pressures,  one  needs  the 
energetics  of  cluster  formation  within  the  bubble  in  the  form  of  an  interaction  potential  function. 
Due  to  the  unavailability  of  estimates  of  many-body  interaction  energies  between  electronically 
excited  stable  clusters  of  He  atoms  with  solvent  He  atoms,  the  present  model  essentially  employs  a 
sum  of  the  two-body  potentials:  it  is  based  on  using  the  diatomic  sets  a,  X  or  c,  X  depending  on 
the  g/u  state  of  He-He*  potential  energy;  thus  it  effectively  ignores  certain  many-body  effects.  It  is 
hoped  that  the  qualitative  description  of  bubble  states  of  high  pressure  liquid  Helium  should  be 
reasonable  with  this  model  although  the  quantitative  nature  of  the  present  results  may  be  changed 
by  such  many-body  energetics. 

Our  Monte  Carlo  (MC)  simulations  are  based  on  the  recently  developed  He2*(a  ^Su+  and 
c  3£g+)  potentials9  in  a  system  consisting  of  one  excited  He*(^S)  atom  situated  in  bulk  liquid  He 
environment,  the  atoms  of  the  latter  interacting  with  each  other  via  the  effective  (labelled  X-state 
here  for  convenience)  Aziz  potential.10 

<p(r)  =  el  be-“R  -  {^R-6  +  Cg  R~8  +  C,^"10} 

exp[-(d/R-l)2],R<d 

R>d  (205) 


f(R)]  (204) 


95 


where  R  =  r/rm,  rm  =  2.9673A,  b  =  544850.4,  a  =  13.353384,  C6  =  L  3732412, 
Cg  =  0.4253785,  =  0. 1781,  d=  1.241314,  e/kg=10.8K,  and  kB  is  the 

Boltzmann  constant. 

The  above  potential  parameters  were  obtained10  by  fitting  to  the  He2  dimer  electronic  and 
vibrational  energy  levels,  second-virial  coefficients,  thermal-conductivity,  high-temperature 
viscosity,  and  differential  cross-section  measurements.  The  Aziz  pair  potential  has  recently  been 
used  in  a  variety  of  condensed  phase  simulation  studies11  and  found  to  reproduce  adequately  a 
sizable  number  of  experimental  data  particularly  the  equation  of  state12,  appropriate  to  our 
moderately  high  pressure  studies13. 

Jordan,  Siddiqui  and  Siska9  have  recently  developed  interaction  potential  models  for  the  a 

-*2  -f- 

J  E  u  and  c  3  £  states  of  He2  on  the  basis  of  generalized  Morse-Hulburt-Hirschfelder  functions 

to  describe  both  the  well  regions  (Region  A)  and  barriers  (Region  B)  and  of  improved  Tang- 
Toennies  form  for  the  long-range  van  der  Waals  interaction  (Region  C).  We  refer  the  reader  to 
their  work  for  the  functional  forms  in  each  region.  The  17-parameter  (see  table  II  of  ref.  4a  except 
that  due  to  a  misprint,  the  sign  of  the  c  parameter  entering  Region  B  interaction  must  be  changed 
for  the  a-state  potential  fit)  models  have  been  obtained  by  using  the  data  from  crossed  beam 
scattering  experiments,  ab  initio  results  and  low-temperature  exchange  rates.  The  a-state  potential 
has  an  attractive  well  for  separations  of  about  lA,  followed  by  a  high  repulsive  barrier  around 
2.7  A.  The  excited  a-state  potential  has  a  hump  in  the  region  of  2.7 A  below  which  is  located  the 
bound  He2*(^£u+)  species  having  an  equilibrium  separation  of  close  to  lA.  Similar  behavior  is 
also  found  for  the  c-state  potential  corresponding  to  He2*(^£g+)  the  main  difference  from  the  a- 

state  being  the  larger  barrier.  The  minima  for  the  c-  and  a-states  are  very  close.  Figure  15  shows  a 
comparison  of  the  intermediate  range  parts  of  the  *,  i^-Pe  interactions  as  well  as  the  Siska  a-  and  c- 
states  which  control  the  nature  of  the  Helium  bubble  u,  -  surrounds  the  excited  atomic  He*  species. 
In  order  to  test  the  effect  of  interaction  potential  on  the  structure  of  the  formed  bubble,  the  available 
new  ab  initio  potential  for  the  long  range  part  of  the  a-state  devised  by  Konowalow  and 
Lengsfield14  (denoted  here  as  the  KL  potential)  has  also  been  employed  in  one  set  of  simulation 
runs13  (see  sec.  IV).  The  well  depth  of  the  KL  potential  is  almost  a  factor  of  4  larger  than  the 
small  van  der  Waals  well  of  the  Siska  a-state  and  the  KL  hump  is  slightly  higher  than  the  Siska 
hump.  Except  for  this  test13  the  potentials  of  ref.  9  were  exclusively  employed  in  this  work. 

3.  Simulation  Details 

The  MC  simulations  were  done  on  a  He*/He  solution  employing  the  above  described 
interaction  potentials  in  a  system  consisting  of  one  excited  He*(-3S)  atom  situated  in  a  bath  of 
manyHe-atoms  with  periodic  boundary  conditions  representing  bulk  liquid  He.  These  simulations 
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basically  pose  two  problems:  (1)  convergence  of  the  results  because  of  the  infinite  dilution  studies, 
and  (2)  trapping  of  He  atoms  (due  to  the  absorbing-sphere  nature  of  the  a-  and  c-state  potentials 
that  allow  a  barrier  crossing  to  occur  resulting  in  a  temporary  absorption  of  solvent  atoms)  and 
formation  of  compact  clusters  of  He  atoms  (that  result  upon  several  solvent  atoms  crossing  over  to 
the  region  below  the  barrier,  that  is  to  the  left  of  the  *  marked  in  Fig.  15).  The  former  problem  has 
been  addressed  by  performing  very  long  and  large  scale  simulations  for  better  statistics.  The 
number  of  particles  needed  for  production  runs  depends  on  the  number  of  solvation  shells  to  be 
included  in  the  simulation  box.  In  order  to  study  the  structure  of  the  second  solvation  shell  at  high 
pressures  (which  is  not  expected  to  be  completely  determined  by  the  256  particle  simulations)  as 
well  as  to  eliminate  simulation  box-edge  effects  on  the  structure  of  the  bubble,  we  performed  500 
particle  simulations  where  necessary.  The  bubble  structure  has  been  obtained  from  an  average  of 
three  MC  runs  at  each  pressure  with  appropriate  MC  move  step-size;  each  run  consisted  of  one 
million  MC  moves.  We  then  adapted  a  nine-point  moving  polynomial  smoother  algorithm15  on  the 
mean  radial  distribution  function  so  as  to  eliminate  numerical  noise  which  can  obscure  relevant 
features  in  the  distribution. 

The  solvent  atom's  access,  on  passing  over  the  a-  or  c-state  barrier,  to  the  attractive  region 

below  the  barrier  adds  to  the  complexity  of  simulations  due  to  trapping  of  solvent  He  atoms  and 
* 

formation  of  an  Hen  complex  in  the  inner  region  of  the  bubble.  This  may  cause  simulation 

problems  not  only  in  the  formation  of  the  bubble  but  also  in  characterizing  the  realistic  structural 

features  of  the  pertinent  bubble  states.  Further,  it  is  expected  that  dynamical  simulations,  necessary 

to  obtain  reliable  GLE  parameters,  cannot  be  easily  carried  out  with  such  barrier-crossing  (see  Fig. 

15  for  comparison  of  energetics)  and  reaction  problems  involving  the  solvent  In  order  to  examine 

these  effects  on  the  structure,  we  concentrated  on  MC  simulations.  Towards  the  goal  of  assessing 

these  barrier  crossing  effects  on  the  bubble  state,  we  experimented  with  different  ways  of 

* 

eliminating  the  trapping  of  He  atoms  to  form  the  Hen  species.  As  a  result,  the  problem  of 

characterization  of  the  He*  bubble  in  liquid  He  is  examined  here  based  on  two  types  (see  below)  of 
MC  simulations  of  an  excited  He*  species  ( a  and  c  states )  in  the  condensed  phase  He  environment 
(solution). 

A.  Straight  Simulation: 

In  this  study,  we  mostly  performed  a  standard  or  natural  (straight)  simulation  of  the 
solution  system  allowing  the  uninterrupted  formation  of  Hen*  complexes.  This  was  done  to 
explore  the  effect  of  a  nonequilibrium  reactive  flux  being  allowed  to  naturally  influence  the 
solvation  structure.  The  simulation  procedure  used  was  standard  except  that  He*  is  placed  at  the 
centre  of  the  simulation  cell  in  the  beginning  of  the  run. 
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Figure  15.  Pairwise-additive  potentials,  U(r),  for  the  ground  (X)  state  of  He2,  and  for  the 
excited  a-  and  c-states  of  He2*.  The  inset  figure  shows  the  r-dependence  of  the 
potentials  for  small  values  of  r.. 
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B.  Prevention  of  clustering: 

Simulation  of  solvent  structure  around  a  specified  metastable  species,  Hen*  requires  that 

the  species  be  somehow  kept  from  transforming  (by  barrier  crossing  reaction  with  the  solvent)  to 
other  possible  species  during  simulation.  We  have  explored  the  effects  of  the  following  two  simple 
schemes  of  preventing  such  a  clustering. 

In  the  simulations  described  below,  we  simplify  the  algorithm  by  fixing  the  He*  atom  at  the 
centre  of  our  simulation  cell  and  keeping  it  fixed  there  (clamped)  during  the  MC  run.  This  was 
done  after  testing  such  a  procedure  in  comparison  to  an  unclamped  one.  Clamping  also  helps 
economize  prevention  of  the  possible  disruption  of  the  formed  bubble,  the  possible  disruption 
being  due  to  the  He*  moving  to  the  edge  of  the  simulation  cell  during  the  simulation. 

/.  Simple  Rejection: 

In  order  to  restrict  the  trapping  of  He  atoms  below  the  barrier  region  by  the  simple  rejection 
method,  we  opted  to  reject  any  attempted  MC  move  if  the  distance  between  the  He*  and  He  is  less 
than  that  of  the  corresponding  barrier  location  (5.1  and  3.7  bohrs  for  a-  and  c-  states  respectively). 
The  number  of  such  rejected  moves  is  monitored  along  with  various  thermodynamic  quantities 
such  as  energy,  and  pressure.  In  order  to  assess  the  effect  of  such  rejected  moves,  a  series  of  MC 
simulations  were  performed  for  both  a-  and  c-states  at  two  pressures.  All  the  simulation  results  are 
obtained  by  averaging  over  2  to  3  million  moves  for  256  particles.  Where  necessary,  we  have 
performed  500  particle  simulations  as  are  mentioned  below. 

2.  Flow  Instead  of  Rejection: 

Though  the  MC  moves  are  being  rejected  in  the  above  simple  rejection  procedure  if  the 
atoms  go  inside  the  barrier  location,  the  solvent  atoms  being  left  at  or  near  the  barrier  location  is 
expected  to  contribute  some  additional  features  in  the  structure.  In  order  to  eliminate  such  effects 
due  to  the  simple  rejection  of  reaction  around  the  barrier  location,  in  the  flow  type  of  solvation 
model,  atoms  trying  to  enter  the  bubble  (i.e.,  either  dimerize  or  form  higher  clusters)  are  instead 
displaced  across  to  the  respective  octant  comers  of  the  cubic  simulation  cell.  It  can  be  visualized 
that  this  results  in  a  flow  system. 

4.  Results  &  Discussion 

To  date,  experimental  studies  on  helium  metastables  and  the  theoretical  models  employed 
for  them  have  been  in  the  low  pressure  and  temperature  regime.  The  occurrence  of  Helium  atoms 
and  He2*  species  have  been  documented  in  this  regime.  Higher  order  clusters  of  He  atoms  such  as 

might  be  expected  at  higher  pressures  are  being  explored  here  for  the  first  time.  The  studies 
peformed  by  us  so  far  are  largely  in  the  nature  of  an  exploration  of  the  technology  required  to 


100 


characterize  these  highly  excited  species  in  condensed  phases.  We  summarize  the  salient  points 
that  have  arisen  in  the  course  of  our  computer  experiments  while  presenting  a  discussion  of  the 
relevant  structural  and  dynamical  issues  many  of  which  are  accessible  to  observation  by  modem 
experiments  not  yet  performed  for  the  present  systems. 

A.  Bubble  Structure 

The  existence  of  a  bubble-like  structure  in  solution  is  clearly  seen  upon  examining  the  radial 
distribution  function  (RDF)  for  the  solute-solvent  pairs  (computed  by  averaging  over  the  MC  runs 
made).  We  discuss  differences  among  the  RDF’s  in  terms  of  the  potential  function  employed  for 
the  solute- solvent  interaction  (a-  or  c-state)  and  other  parameters  of  the  simulation  (eg.  pressure, 
imposed  boundary  conditions  (cf:  sec.  Ill)  etc.). 

To  begin  with,  we  tested  the  effect  of  clamping  the  He*  atom  on  the  RDF's.  Essentially 
identical  results  were  obtained  for  the  c-state  bubbles  at  pressures  of  0.5  and  1.4  GPa  by 
constraining  or  not  constraining  the  He*  (but,  in  the  latter  case,  ensuring  that  the  box  is  large 
enough  for  the  entire  bubble  to  stay  inside  it  for  the  whole  simulation  run).  Secondly,  The  KL  a- 
state  bubble  was  simulated  to  see  if  the  differences  in  the  two  a-state  potentials  (of  ref.  9  and  14) 
are  detectable  in  the  bubble  structure  and  hence  judge  the  adequacy  of  the  long-range  potential 
being  employed.  Although  minor  quantitative  differences  in  peak  heights  were  noticeable,13  the 
peak  positions  were  found  to  be  essentially  the  same.  These  changes  in  peak  heights  were 
qualitatively  consistent  with  the  difference  in  the  potentials  employed:  these  are  a  shorter  first  peak 
attributable  to  the  higher  barrier  of  KL  and  a  minor  shift  of  probability  to  larger  distances  (closer  to 
the  van  der  Waals  minimum),  but  not  significant  accumulation  to  reflect  any  important  role  due  to 
the  deeper  well  there.  Thus  the  a-state  potential  of  Ref.  9  was  retained  for  the  rest  of  the  work. 

Figure  16  presents  our  results  for  the  He*-He  radial  distribution  functions  (RDF)  of  c-state 
at  0.5,  1.4,  and  10  GPa  obtained  using  the  Straight  model  simulations.  The  RDF's  for  0.5  and 
1.4  GPa  show  no  density  of  He  atoms  below  the  peak  structure  located  at  5-63Q.  It  is  seen  by 

comparing  to  the  distance  scale  of  the  potential  functions  in  Fig.  15  that  the  structured  RDF  for 
these  two  pressures  represent  the  density  distribution  of  the  first  few  solvation  shells  of  the 

iji 

"bubble"  within  which  the  atomic  He  species  is  located.  The  large  RDF  peak  for  10  GPa  near 
2ag  can  be  understood  to  reveal  the  formation  of  a  possible  cluster,  Hen  j  and  the  remaining  outer 

structure  as  due  to  the  formation  of  the  bubble.  The  clear  demarkation  between  the  two  (inner 
versus  outer  peaks)  seen  in  the  fact  that  the  probability  becomes  zero  near  about  4  bohrs  and  rises 
slowly  from  then  onwards  indicates  that  it  makes  sense  to  call  the  object  responsible  for  the  short 
range  structure  a  "bubble  state"  just  as  was  done  for  the  atomic  and  dimer  cases  at  lower  pressures 
based  on  experiments.  It  may  be  asked  if  the  structure  of  the  bubble  solvation  shells  will  be 
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Figure  16.  The  structure  of  the  He*(c)  bubble  as  a  function  of  pressure: 

_ 10.0  GPa 

_  1.4  GPa 

_  0.5  GPa 
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influenced  due  to  the  additional  interactions  of  the  solvent  with  the  absorbed  solvent  atoms  present 
in  the  cluster  species  within  the  bubble.  If  these  interactions  are  important,  the  solvation  structure 
at  these  pressures  would  correspond  to  a  different  bubble  state  for  different  cluster  species.  We 
have  found  that  (see  below)  in  two  of  the  lower  pressures  studied  here  (0.5  and  1.4  GPa),  this  is 
probably  not  a  significant  effect  although  this  effect  could  be  present  at  10  GPa.  The  average 
number,  <n>,  of  atoms  present  in  the  bubble  have  been  monitored  during  the  simulations  and  are 
included  in  Table  4. 

Figure  17  contains,  for  the  same  pressures,  the  corresponding  RDFs  from  simulations 
based  on  the  a-state  interaction.  Comparison  of  Figs.  17  and  18  indicate  the  following:  (1)  there  is 
a  shift  in  peak  positions  for  both  a-  and  c-states  towards  smaller  radii  for  higher  pressures.  The 
shift  in  the  peak  positions  as  well  as  the  appearance  of  a  second  solvation  shell  at  higher  pressure 
can  be  readily  understood  in  terms  of  packing  effects  based  on  available  volume;  (2)  the  bubble 
peak  positions  do  not  depend  on  the  g/u  electronic  state  characters  of  the  He*-He  interactions  and 
could  be  because  they  basically  lie  close  in  energy;  (3)  the  difference  in  peak  heights  (higher  for  the 
c-state  bubble  at  10  GPa)  could  be  due  to  a  combination  of  the  higher  slope  of  the  c-state  barrier 
near  5.7a0  as  well  as  excluded  volume  forces  from  additional  atoms  (<n>=3.6  for  c-state  in  Table 

4)  of  the  cluster  within  the  bubble.  Clearly,  at  10  GPa,  the  a-state  has  more  particles  within  the 
bubble  (see  <n>  =  6.1  for  a-state  in  Table  4)  at  this  pressure  but  they  are  found  unable  to 
contribute  enough  additional  repulsion  to  compensate  also  for  the  higher  slope  of  c-state  potential 
near  5.7aQ;  and  (4)  There  is  significant  cluster  formation,  revealed  by  short-range  structures  in  the 

a-bubble  RDF,  at  the  two  lower  pressures,  in  contrast  to  the  c-bubble.  The  fact  that  the  outer  part 
of  the  RDF's  for  the  two  states  at  these  lower  pressures  (0.5  and  1.4  GPa)  coincide,  just  as 
expected  from  the  almost  degenerate  nature  of  the  g/u  potentials  for  distances  beyond  6slq,  confirms 

that  the  potential  energy  contribution  of  the  additional  atoms  from  the  cluster  located  within  the 
bubble  to  solvation  structure  modification  is  negligible  at  least  when  <n>  is  as  high  as  around  6 
(Note  that  although  the  <n>  value  for  c-state  cluster  at  10  GPa  in  relation  to  point  3  above  is  only 
between  3  and  4,  the  peak  position  in  question  is  of  shorter  range  at  10  GPa). 

B.  Clustering  Propensity 

The  value  of  <n>  is  seen  to  depend  on  which  electronic  state  is  involved  so  that  electronic 
nonadiabaticies  would  have  an  effect  on  this  simple  quantity.  The  propensity  for  the  formation  of 
Hen*  complexes  was  found  more  pronounced  in  the  a-state  compared  to  the  c-state  and  is 
attributable  to  their  barrier  height  difference.  In  this  work  we  found  high  pressure  values  of  <n> 
(see  Table  4)  significantly  larger  than  the  value  of  2  that  presumably  occurs  at  ordinary  low 

pressures  (as  in  reports  based  on  previous  experiments).  Table  4  also  presents  additional 

<¥ 

simulation  data  for  a  and  c  states,  such  as  the  average  He  -He  bond  separation  in  the  cluster  that  is 
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trapped  within  the  bubble.  A  unit  value  for  RDF  corresponds  to  the  uncorrelated  limit  of  average 
liquid  bulk  pair  probability.  The  pair  probability  for  the  inner  part  (i.e.,  within  the  bubble)  of  the 
kDF  is  found  to  be  an  order  of  magnitude  larger  (with  peak  heights,  seen  for  e.g.,  as  in  the  inset 
of  Fig.  17  shown,  of  around  10-30)  than  the  average  long-range  bulk  pair  probability.  These 
relative  RDF  magnitudes  reflect  the  much  more  compact  nature  of  the  cluster  located  within  the 
bubble  when  compared  to  the  liquid  bulk.  This  compact  nature  leads  us  to  a  physical  picture  akin 
to  gas/surface  type  collision  dynamics  between  the  species  in  the  bubble  and  the  concave  inner 
surface  of  the  bubble. 

The  <n>  values  for  a-bubble  appear  to  have  saturated  to  a  value  around  6  at  the  highest 
pressures  studied  here,  and  an  order  of  magnitude  increase  in  the  pressure  beyond  the  saturation 
point  is  seen  to  maintain  the  bubble  structure  intact.  What  <n>  reveals  is  the  relative  change  in  the 
free  energy  of  cluster  formation  within  the  bubble.  The  variation  of  <n>  may  be  termed  a  static 
solvent  effect.  It  must  be  noted  here  that  MC  simulations  do  not  yield  the  n-value  dynamical 
fluctuations  which  may  affect  the  observability  conditions  for  the  n-mers. 

The  questions  of  how  much  these  clustering  reactions  influence  the  bubble  and  the 
dynamics  of  a  given  n-mer  in  solution  are  not  addressed  by  these  MC  simulations  within  the 
Straight  model.  In  attempting  to  develop  some  insight  into  these  processes,  we  have  developed  the 
other  two  models  described  in  Subection  3  that  prevent  further  clustering  of  a  specified  n-mer  with 
a  solvent  atom  to  form  an  (n+l)-mer.  Results  for  RDFs  from  the  Flow  model,  the  Simple 
Rejection  model,  and  the  Straight  model  for  the  a-state  at  1.4  GPa  are  compared  in  Fig.  18.  The 
first  two  types  of  simulations  (  termed  Monomer'l  and  Monomer'2  in  Fig.  18)  constrain  the 

jft  * 

bubble  species  to  be  solely  He  atoms  by  preventing  the  formation  of  even  He2  .  Clearly,  neither 

the  peak  positions  nor  the  peak  structural  features  are  in  good  agreement  with  each  other  in  Fig. 

1 8,  thereby  displaying  tremendous  sensitivity  to  the  simulation  model  choice  and  the  important  role 

of  reaction  in  bubble  structure  maintenance,  at  least  for  this  n=l  case.  This  <n>  is  quite  low 

compared  to  the  average  <n>  =  6.3  value  of  Table  4  for  this  pressure.  Thus  the  atomic  species 

examined  is  probably  an  unstable  one  under  the  conditions  of  1.4  GPa.  Naturally,  no  rejection 

* 

scheme  makes  physical  sense  unless  the  level  n,  at  which  Hen  is  fixed,  is  a  stable  (albeit 

metastable)  species  in  the  condensed  phase  matrix. 

We  conclude  the  discussion  of  the  nonclustering  model  studies  with  a  few  remarks  on  the 
technology  used  to  do  the  simulations.  It  was  found  that  the  256-particle  Flow  model  simulation 
resulted  in  RDF  values  greater  than  1.0  around  the  boundary  of  the  simulation  cell.  In  order  to 
avoid  this  spurious  effect  as  well  as  to  reduce  any  possible  effect  on  the  main  peak  of  the  bubble, 
we  performed  500-particle  simulations  and  this  is  the  one  included  in  Fig.  1 8.  The  first  peak 
position  and  height  are  reproducible  in  both  sets  (256  and  500  particle  runs)  of  simulations.  The 
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Figure  17.  The  structure  of  the  He*(a)  bubble  as  a  function  of  pressure: 

_ 10.0  GPa 

_  1.4  GPa 

_  0.5  GPa 


The  inset  shows  the  short-range  structure. 


1.8 
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Figure  18.  The  structure  of  the  He* (a)  bubble  for  different  simulation  strategies  that  prevent 
clustering  with  the  solvent  (as  discussed  in  Subsections  3  and  4  of  this  section). 
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Flow  model  RDF  in  Fig.  18  may  imply  that  although  the  bubble  may  have  been  formed  it  is 
disrupted  or  diffused  out  because  the  particles  responsible  for  its  formation  are  being  taken  out  to 
the  comers  of  the  cell.  Clearly,  the  simple  prescriptions  of  the  nonclustering  models  employed 
here  are  insufficient  for  achieving  the  purpose  of  realistic  characterization  of  n-mers  when  n  is  very 
different  from  <n>.  It  is  hoped  that  improved  simulation  schemes  can  be  devised,  but  (we  feel) 
perhaps  only  for  situations  when  n  is  closer  to  <n>  so  as  to  study  pertinent  species  (which  should 
not  have  n  that  is  too  different  from  <n>). 

5.  Conclusion 

We  have  presented  results  of  a  thermodynamic  state  dependent  MC  study  of  the  Helium 
bubble  states  arising  from  electronically  excited  He  metastables  being  trapped  in  a  condensed 
phase  He  environment.  Even  in  this  classical  regime,  bubbles  were  found  to  exist  and  house  the 
excited  species,  although  the  bubble  contents  are  a  sensitive  function  of  pressure.  Various  RDF’s 
were  presented  to  analyze  and  interpret  the  bubble  changes  with  pressure,  electronic  state  features 
etc.. 

This  type  of  simulation  is  also  highly  pertinent  to  materials  problems  such  as  embrittlement 
in  fusion  reactor  materials  associated  with  generation  of  He  that  have  attracted  considerable  interest 
in  recent  years.  Bombardment  with  heavy  noble-gas  ions  is  widely  used  in  various  modem 
preparation  processes16  for  metal  surfaces  and  films.  It  has  been  shown  that  the  He  embrittlement 
observed  in  metals  was  due  to  the  trapping  and  subsequent  bubble  formation  at  radiation  induced 
defects  as  well  as  in  the  absence  of  radiation  damage17.  To  evaluate  the  effects  of  noble  gases  on 
physical  or  material  properties  of  metals,  it  is  therefore  important  to  be  able  to  characterize  noble- 
gas  bubbles  and  to  develop  an  understanding  of  their  behavior.  There  have  been  a  number  of 
atomistic  calculations4  to  explain  the  bubble  formation  and  the  involved  self-trapping  mechanism  of 
He  in  metals.  We  have  revealed  certain  technically  feasible  aspects  as  well  as  raised  some  pertinent 
issues  relating  to  characterizing  the  role  of  a  reactive  species  being  present  within  bubbles. 
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Table4a:  MC  Results  of  He*(23S)  in  He  (1JS)  bath  at  300K 


'n-mer'  Structure 


v  m/Cm3  mol"1 

(P)/GPa 

(U)/kJ  mol*1 

(rHe*-lfe)/a0 

<n> 

17.95 

0.29 

0.36 

2.29 

5.6 

13.17 

0.55 

-2.0 

2.29 

5.9 

0.57 

0.85 

— 

0 

9.02 

1.4 

-0.9 

2.32 

6.3 

1.5 

2.2 

— 

0 

4.76 

10.1 

10.4 

2.24 

6.1 

10.3 

11.9 

2.14 

3.6 

4.35 

13.8 

14.2 

2.27 

6.1 

13.9 

15.6 

4.2 

% 

(a)  First  and  second  row  values  for  each  volume  correspond  to  a-and  c-states  of  He-He 
interaction. 


Ill 


REFERENCES 


1  (a)  W.  S.  Dennis,  E.  Durbin,  Jr.,  W.A.  Fitzsimmons,  O.  Heybey,  and  G.  K.  Walters,  Phys. 
Rev.  Letts.  23,  1083  (1969);  (b)  J.  C.  Hill,  O.  Heybey,  and  G.  K.  Walters,  ibid.  26,  1213 
(1971);  (c)  John  L.  Watkins,  Jonas  S.  Zmuidzinas,  and  Gary  A.  Williams  Physica  108B, 
1313  (1981). 

2  (a)  See  (la);  (b)  A.  P.  Hickman  and  N.  F.  Lane,  Phys.  Rev.  Letts.  26,  1216  (1971);  (c)  A. 
P.  Hickman,  W.  Steets,  and  N.  F.  Lane,  Phys.  Rev.  B12,  3705  (1975);  Y.  M.  Shih  and  C- 
W.  Woo  Phys.  Rev.  A8,  1437  (1973). 

3  (a)  E.  V.  Komelsen,  Can.  J.  Phys.  48,  2812  (9170);  (b)  W.  D.  Wilson  and  C.  L.  Bisson 
Radiad.  Eff.  22,  63  (1974). 

4  (a)  W.  D.  Wilson,  C.  L.  Bisson,  and  M.  I.  Baskes  Phys.  Rev.  B  24,  5616  (1981);  (b)  Refs, 
in  Advances  in  the  Mechanics  and  Physics  of  Surfaces,  ed.  R.  M.  Latanision  and  T.  E. 
Fisher,  Harwood  Academic  Publishers,  (1983);  (c)  K.  O.  Jensen  and  R.  M.  Nieminen,  Phys. 
Rev.  B  35,  2087  (1987);  Ibid.  36,  8219  (1987). 

5  J.  P.  Hansen  and  E.  L.  Pollock,  Phys.  Rev.  A5,  2214  (1972). 

6  OJgQ  and  ujc1  are  related  to  the  second  and  fourth  frequency  moments  of  the  system  spectral 

density  and  are  defined  for  eg.  by  S.  A.  Adelman  J.  Chem.  Phys.  73,  3145  (1980);  (b)  For 
definition  of  GLE  parameters  in  a  partial  clamping  approximation  scheme,  see  S.  A.  Adelman 
J.  Chem.  Phys.  81,  2776  (1984). 

7  D.  A.  Micha,  J.  Chem.  Phys.  78,  7138  (1983). 

8  P.  K.  Swaminathan,  B.  C.  Garrett,  and  C.  S.  Murthy  J.  Chem.  Phys.  88,  2822  (1988);  B. 
C.  Garrett,  P.  K.  Swaminathan,  C.  S.  Murthy,  and  M.  J.  Redmon  J.  Chem.  Phys.  87,  3207 
(1987). 

9  (a)  R.  M.  Jordan,  H.  R.  Siddiqui,  and  P.  E.  Siska  J.  Chem.  Phys.  84,  6719  (1986);  (b)  R. 
M.  Jordan  and  P.  E.  Siska  J.  Chem.  Phys.  80,  5027  (1984). 

10  R.  Aziz,  V.  P.  S.  Nain,  J.  S.  Carley,  W.  L.  Taylor,  and  G.  T.  McConville,  J.  Chem.  Phys. 
70,  4330  (1979). 

11  (a)  D.  Levesque,  J.  J.  Weis,  and  M.  L.  Klein,  Phys.  Rev.  Lett.  51,  670  (1983);  (b)  D. 
Levesque,  J.  J.  Weis,  and  P.  Loubeyre,  Phys.  Rev.  B  34,  178  (1986);  (c)  D.  Frenkel,  Phys. 
Rev.  Lett.  56,  858  (1986);  (d)  D.  M.  Ceperly  and  H.  Partridge,  J.  Chem.  Phys.  84,  820 
(1986);  (e)  P.  Loubeyre,  Physica  139&140  B,  224  (1986);  (f)  P.  Loubeyre,  Phys.  Rev. 
Lett.  58,  1857  (1987);  (g)  M.  Takahashi,  J.  Phys.  Soc.  Japan,  55,  1952  (1986);  (h)  W.  B. 
Smith  and  K.  Singer  Molec.  Phys.  (in  press). 

12  R.  L.  Mills,  D.  H.  Liebemnberg,  and  J.  C.  Bronson,  Phys.  Rev.  B  21,  5137  (1980). 

13  P.  K.  Swaminathan  and  C.  S.  Murthy  (unpublished  results). 

14  D.  Konowalow  and  B.  Lengsfield,  J.  Chem.  Phys.  87,  4000  (1987). 


112 


15  A.  Savitsky  and  M.  J.  E.  Golay,  Anal.  Chem.  36,  1627  (1964). 

16  (a)  R.  G.  Musket,  W.  McLean,  C.  A.  Colmenares,  D.  M.  Makowiecki,  and  W.  J.  Siekhaus, 
Appl.  Surf.  Sci.  10,  143  (1982);  (b)  S.  Matteson  and  M.  A.  Nicolet,  Annu.  Rev.  Mater.  Sci. 
13,  339  (1983). 

17  (a)  G.  J.  Thomas,  W.  A.  Swansiger,  and  M.  I.  Baskes,  J.  Appl.  Phys.  50,  1942  (1979);  (b) 
G.  J.  Thomas  and  R.  Bastasz,  J.  Appl.  Phys.  52,  6426  (1981). 


113 


Heatbath  Models  for  Helium  Bubble  States 

L  Introduction 

Production  of  4  He,  the  identification  of  (metastable)  excited  states  of  4  He  atoms,  and 
exploration  of  associated  problems  in  fusion  reactors1  have  all  been  active  areas  of  experimental 
study  for  some  time.  On  the  theoretical  front,  in  the  use  of  ab  initio,  lattice  dynamical,  and 
statistical  mechanical  approaches,  helium  has  long  attracted  many  researchers  because  it  constitutes 
a  fundamental  system  that  is  challenging  for  technique  developments  involving  quantum  effects  in 
various  contexts.2 

The  indications,  obtained  in  recent  high-pressure  experiments3  using  the  diamond-anvil¬ 
cell  technique,  of  a  possible  new  solid  phase  around  300  K  along  the  melting  curve,  have  resulted 
in  extensive  theoretical  research4  effort  on  solid  He.  This  is  largely  because  of  the  availability  of  a 
reasonably  accurate  pair  potential20  and  the  recent  advances  in  the  computer  simulation  techniques 
for  studying  phase  transitions5  and  estimating  ffee-energies.6  There  have  also  been  a  few  studies 
in  the  fluid  phase,7  dealing  with  the  thermodynamic  properties  and  equilibrium  structure. 

Our  goal  is  to  theoretically  treat  the  dynamics  of  electronically  excited  species  and  the 
condensed  phase  influence  on  such  species.  We  have  chosen  4He  for  the  benchmark  studies 
necessary  for  computational  method  developments.  In  the  present  work,  we  performed  classical 
Monte  Carlo  (MC)  and  classical  molecular  dynamics  (MD)  (where  appropriate)  simulations  on 
liquid  He  as  well  as  solution  systems  consisting  of  one  excited  He*(3S)  atom  situated  in  bulk 
liquid  He  environment  Our  work  is  focused  on  state  points  near  300  K  in  the  high  pressure  (GPa) 
regime.  We  have  not  addressed  the  quantum  translational  aspect  and  the  three-body  force-effects 
on  the  structure  and  dynamics  of  liquid  He,  although  they  deserve  a  special  and  large  effort  in 
order  to  accurately  quantify  such  effects.  Although  these  effects  are  at  the  outset  expected  to  be 
important  for  He,  it  is  anticipated  from  the  existing  investigations  on  He8  that  quantum  effects  have 
only  a  minor  role  to  play  at  high  densities  (high-pressure  10-20  GPa).  Our  work  here  is  based  on 
classical  simulations  using  the  recently  developed  Aziz  potential.20 

We  now  turn  to  the  question  of  the  degree  of  suitability  of  the  Aziz  effective  pair  potential 
for  realistic  simulations  of  4He  and  effects  of  three-body  forces.  Although  the  Aziz  pair  potential 
has  been  widely  used  in  the  recent  high  density  studies  on  4He  which  are  referred  to  above,  there 
remain  only  two  recorded  deficiencies9  of  the  model:  it  is  (1)  too  repulsive  for  distances  less  than 
1.8  A  and  (2)  fails  to  explain  the  thermodynamic  stability  of  the  bcc  phase  around  300  K.  The 
former  weakness  of  the  model  will  matter  only  at  very  high  densities  (pressures  above  60  GPa). 
The  latter  may  well  be  due  to  the  neglect  of  three-body  forces.  It  has  in  fact  been  recently  shown10 
that  self-consistent  phonon  and  MC  calculations,  involving  Aziz's  pair  potential  with  three-body 
interactions,  produce  results  in  excellent  agreement  with  the  experimental  equation  of  state  at  high 
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pressures.1 1  It  should  however  be  noted  that  such  refinements  are  only  necessary  for  resolving  the 
stability  of  equally  competing  solid  phases  such  as  fee  and  hep  and  perhaps  for  detailed  description 
of  dynamical  properties. 

The  present  simulations  of  liquid  solvation  structure  can  yield  the  associated  heat-bath 
modelling  MTGLE  parameters  for  realistic  quenching  studies  (e.g.,  in  He-He*  matrix  studies). 
This  involves  certain  approximate  statistical  mechanical  relations  between  structural  properties  and 
the  frequency  moments  of  the  GLE  heatbath  spectral  density.  Such  prescriptions  for  determination 
of  parameters  for  stochastic  dynamical  treatments  are  delineated  in  a  series  of  papers  by  Adelman 
and  coworkers,12  in  particular,  within  the  partial  clamping  formalism  developed  for  solvated 
polyatomic  entities.  A  dynamical  simulation  of  solvation  is  essential  for  the  exact  characterization 
of  the  dynamics  of  a  heatbath  since  the  real  time-dependent  friction  kernels  are  dynamical 
quantities.  In  this  report,  following  Adelman  and  coworkers,  we  examine  the  validity  of  a 
gaussian  model  for  the  time  dependent  solvent  friction,  and  having  obtained  some  MD-based 
information  for  the  degree  of  validity  of  the  gaussian  model  for  the  solvent,  we  apply  the  gaussian 
model  to  the  MC-based  structural  (radial  distribution  functions)  data  from  a  solution  of  He*  in  He 
(reported  earlier).13  This  kind  of  approach  is  necessary  especially  because  heatbath 
parameterization  via  MD  simulations  is  difficult  for  the  solution  problem  due  to  the  possibility  of  a 
reaction  between  He*  and  the  solvent  He.  From  our  present  calculations,  we  obtain  evidence  for  a 
decrease  in  pressure  sensitivity  upon  electronic  excitation  in  the  high  pressure  He-matrix,  for  a 
modelled  excited  He*  state  (specifically,  we  employed  the  c-state  for  He-He*  interaction) 
compared  to  ground  state  He  (leading  to  the  Aziz  potential  of  He-He  interaction). 

Subsection  2  describes  the  adopted  interaction  potential  models  and  simulation  details  along 
with  a  brief  description  of  how  the  MTGLE  parameters  have  been  estimated  in  the  present  work. 
Subsection  3  contains  the  main  results  and  discussion. 

A.  Potential  Models  and  Simulations 

We  employed  the  recently  developed  multiparameter  Aziz  model2c  of  effective  pair 
interactions  in  condensed  phase  Helium  in  our  studies  of  the  structure  and  dynamics  of  liquid  He. 
This  effective  pairwise-additive  potential  has  recently  been  employed  in  a  variety  of  condensed 
phase  simulation  studies,4  as  discussed  in  the  introduction.  It  not  only  reproduces  adequately  a 
sizable  number  of  experimental  data  but  also  we  find  that  it  accounts  particularly  well  for  the 
equation  of  state11  appropriate  to  our  present  moderately  high  pressure  studies.  The  functional 
form  of  the  model  used  can  be  written  as 

<P(r)  =  e[  »c-aR-{c6R-6  +  C8R-8  +  C10R-,°}  f(R)) 
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where  R=r/rm,  rm  =  2.9673  A,  b  =  544850  .4,  a  =  13.353384,  C6=  1.3732412, 

C8=  0.4253785,  C1Q=  0.1781,  d  =  1.241314,  e/kB=10.8K,  and  kg  is  the  Boltzmann 

constant.  The  corresponding  potential  parameters  were  determined20  by  fitting  to  the  differential 
cross-section  measurements,  second-virial  coefficients  of  the  gaseous  He,  thermal-conductivity, 
high-temperature  viscosity,  and  He2  dimer  electronic  and  vibrational  energy  levels. 

The  solution  simulations  involve  a  system  consisting  of  one  excited  He*  atom  situated  in 
bulk  liquid  He  environment,  the  atoms  of  the  latter  interacting  with  each  other  via  the  effective  Aziz 
potential,  whereas,  for  the  He* -He  interactions,  we  adopted  the  recently  developed  He2*(a  3£u+ 

and  c  3Ig+)  potentials  of  Jordan,  Siddiqui  and  Siska.14  They  employed  generalized  Morse- 
Hulburt-Hirschfelder  functions  to  describe  both  the  well  regions  (Region  A)  and  barriers  (Region 
B)  and  Tang-Toennies  form  for  the  long-range  van  der  Waals  interaction  (Region  C).  We  refer  the 
reader  to  their  papers  for  the  detailed  functional  forms  in  each  region.  The  17-parameter15  models 
have  been  obtained  by  using  the  data  from  crossed  beam  scattering  experiments,  ab  initio  results 
and  low-temperature  exchange  rates. 

B.  MTGLE  Parameters 

The  MTGLE  parameters  needed  for  the  stochastic  dynamical  modelling  in  liquid  state 
include  the  characteristic  Einstein  frequency,  denoted  co^,  the  coupling  constant,  denoted  cocl) 

and  the  friction  kemal,  denoted  (3(t).  co^  enters  the  stochastic  equations  of  motion  and  governs  the 

initial  energy  transfer  between  the  solute  species  and  (approximately)  the  first  solvation  shell 
species  of  the  solvent  degrees  of  freedom.  coeQ  determines  the  encounter  frequency  for  solute- 

solvent  collisions  that  are  responsible  for  maintaining  or  quenching  the  metastable  species  trapped 
in  bulk.  toCl  controls  the  overall  efficiency  of  energy  transfer  from  the  reagents  to  the  heatbath. 

The  detailed  theory  of  MTGLE  and  expressions  for  the  basic  MTGLE  parameters,  and 

coCl,  in  various  regimes  (exact  as  well  as  for  both  isolated  and  interacting  solute  atoms)  are 

discussed  at  length  by  Adelman  and  coworkers.12  For  completeness,  here  we  give  the  appropriate 
(isotropic  liquid  satisfying  virial  theorem)  formulae  used: 


(208) 


For  pairwise-additive  potentials  between  the  atom- solvent  species,  00  e  and  C0c  ^ ,  in  the 

linear  response  approximation,  can  be  expressed12  in  terms  of  the  corresponding  radial  distribution 
function  (RDF),  g(r): 
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w4c  t  ^  Tr  2g(r)  [{(p" (r)}2  +  {  2cp*  (r)/  r}2]  dr  (210) 

where  m  and  M  are  the  masses  of  the  solute  and  solvent  species,  p  is  the  solvent  number  density, 
g(r)  is  the  solute-solvent  radial  distribution  function,  cp(r)  is  the  interaction  potential  between  the 
solute  and  solvent  species. 

The  time-dependent  friction,  (3(t),  determines  the  frictional  drag  on  the  various  different 
frequencies  participating  in  the  dynamics.  If  one  examines  the  frequency  dependent  function  p(co), 
obtained  by  fourier  transforming  P(t),  the  values  of  P(co)  at  the  frequencies  corresponding  to 
critical  dynamical  regions  on  the  potentials  play  a  key  role  in  describing  important  solvent  effects. 
It  is  possible  to  compute  P(t)  from  the  present  MD  runs  by  employing  the  Volterra  integral 
equation16  relating  various  time  correlation  functions  obtainable  from  the  MD  simulation 
trajectories.  In  addition  to  the  VACF,  these  correlation  functions  include  position  autocorrelation 
as  well  as  position- velocity  cross-correlation  functions.17 

3.  Results  and  Discussion 

Certain  features  of  condensed  phase  Helium  are  unique  and  noteworthy  in  the  simulation 
results.  First  of  all,  the  classical  MC/MD  simulation  procedures  yield  reasonably  converged  results 
for  the  liquid  state  condensed  Helium  only  when  the  number  of  particles  in  the  system  is  near  256 
for  lower  pressures,  and  500  for  the  higher  pressures  (>10  GPa).  This  is  because  there  are  three 
significant  peaks  in  the  radial  distribution  function  (see  below)  and  that  means  two  solvation  shells 
around  each  He  atom,  which  requires  the  larger  number  of  atoms  than  the  108  particles  usually 
employed  for  the  Argon  system.  Further,  the  behavior  of  He  in  the  condensed  phase  is  quite 
different  from  the  case  of  Argon.  The  most  striking  aspect  of  this  is  the  predominant  role  of 
entropy  effects  in  the  liquid  free  energy  of  binding  (this  is  reflected  in  the  fact  that  the  internal 
energy  of  binding  is  an  excess  positive  quantity  as  actually  seen  in  the  simulations  as  well  as  in  the 
configurational  integral,  which  may  be  readily  computed  from  the  known  radial  distribution 
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functions).  In  Table  5,  we  report  the  molar  volumes  at  which  simulations  have  been  performed, 
the  configurational  internal  energies,  observed  pressures,  and  the  factor  PV/RT. 

The  computed  trends  in  liquid  structure,  dynamics,  and  the  MTGLE  heatbath  parameters 
for  pure  liquid  Helium  as  a  function  of  pressure  are  presented  in  Figs.  19  through  21  respectively. 
Liquid  Helium  may  be  treated  classically  in  the  examples  shown,  all  of  which  correspond  to  a 
temperature  of  around  300  K.  As  pressure  is  increased  quantum  effects  become  less,  and  may 
possibly  become  relatively  important  for  the  present  pressures  of  0.3, 0.5  and  1.4  GPa.  This  may 
be  borne  out  from  the  factor  PV/RT,  reported  in  the  last  column  of  table  5.  On  the  basis  of  an 
oscillator  model  obeying  the  Gruneisen  equation  of  state  for  the  solid  and  using  estimated 
Gruneisen  parameter  and  the  Debye  temperature  at  room  temperature  and  high  pressures,  Levesque 

et  al,4a  estimated  AP  qV  /  RT  =  L,  where  AP  ^  is  the  leading  quantum  contribution  to  the 

pressure.  As  our  temperature  and  pressure  range  is  analogous  to  their  studies  along  the  melting 
curve,  one  may  assume  a  value  similar  to  theirs  for  the  quantum  contribution  to  the  pressure.  Such 
a  value  for  the  quantum  contribution  to  the  pressure  may  be  compared  with  the  PV/RT  values 
reported  in  table  5. 

Figure  19  displays  the  RDFs  for  six  different  pressures,  three  of  which  (0.3,  0.5  and  1.4 
GPa)  are  close  to  the  gaseous  fluid  and  three  (10.4,  12.7,  and  14.4  GPa)  are  close  to  a  liquid 
under  high  pressure.  The  RDFs  reflect  the  decreasing  nearest  neighbor  distances  with  increasing 
pressure  and  the  more  diffuse  distributions  of  neighbors  in  the  gaseous  fluid  phase.  The  lOGPa 
and  14  GPa  results  were  obtained  from  500  particle  simulations  and  are  more  accurate  than  the  256 
particle  12.7  GPa  result.  The  He  atom's  velocity  autocorrelation  functions  (VACF)  shown  in  Fig. 
20  display  a  dramatic  change  in  dynamical  character  as  a  function  of  pressure.  In  fact,  the  low 
pressure  VACFs  indicate  gas-like  behavior  with  a  long  hydrodynamic  positive  tail  whereas  the 
high  pressure  ones  show  the  oscillatory  str-  ture  characterizing  a  liquid-like  condensed  phase 
correlation. 18 

Although  the  quantum  effects  are  expected  to  contribute  to  the  structure  and  dynamics  of 
liquid  4He  at  0.3, 0.5  and  1.4  GPa  on  the  basis  of  the  discussions  above,  recent  semiclassical  and 
path-integral  simulation  studies19  on  modelled  rare  gas  atoms  indicate  that  the  effect  of 
semiclassical/quantum  corrections  to  the  RDF  is  not  significant.  The  first  peak  of  RDF  gets 
slightly  broadened  resulting  in  reduction  of  the  peak  height  and  closer  contact  radius  value.  Such  a 
detailed  comparison  for  the  VACFs  from  classical  vs  semiclassical/quantum  simulations  has  not 
been  made  to  our  knowledge.  The  self  diffusion  constants  obtained  from  the  area  of  the  VACFs 
given  by  semiclassical  (Gaussian  wave  packet  simulations)18®  calculations  are  lower  than  are 
obtained  by  classical  MD  simulations.  Thus  further  work  is  needed  to  establish  the  quantum 
effects  on  the  initial  decay  and  the  behavior  of  VACF  at  intermediate  times.  We  can  then  assess  the 
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corresponding  effects  on  MTGLE  parameters.  As  the  quantum  contributions  have  modest 
influence  on  the  RDF  as  discussed  above,  we  may  expect  that  the  quantum  effects  on  co^  would  be 

less  severe  than  on  coCl(see  below)  and  higher  moments. 

The  MTGLE  parameters  have  been  obtained  for  the  tagged  solute  He  in  liquid  He  both  by 
direct  ensemble  average  during  the  simulation  run  itself  using  Eqs.  (206)  and  (207)  (at  the  expense 
of  more  CPU  time  for  the  additional  computation  in  the  forces  loop)  and  from  the  cp(r)  and  g(r)  in 
Eqs.  3  and  4.  For  a  pressure  of  0.3  GPa,  the  a>eQ  and  outvalues  are  18  and  28  THz  from 

ensemble  averages  compared  to  16  and  33  THz  estimated  from  Eqs.  (208)  and  (209).  Similar 
comparison  can  be  seen  from  the  results  in  table  6.  Figure  21  shows  values  for  the  moments20  of 
the  spectral  density,  coeQ  and  coCl,  as  a  function  of  pressure  obtained  as  ensemble  averages.  It  is 

seen  that  these  bulk  liquid  values  are  quite  sensitive  to  (and  exponentially  dependent  on)  the 
pressure.  The  sensitivity  of  the  dynamics  shown  by  the  velocity  autocorrelation  functions  in  Fig. 
20  reflect  the  basis  for  the  MTGLE  parameter  variation^  ell. 

Table  6  gives  the  values  of  the  MTGLE  parameters  coer)  and  <x>Cl  to  be  used  if  the  He*  atom 

is  treated  as  if  it  were  a  diffusing  particle  in  the  liquid  state  with  a-  and  c-  state  solute-solvent 
interactions  assumed.  The  table  also  contain.;  for  comparison  the  corresponding  parameters  for  the 
ground  level  He  atoms  (with  X-state  interaction)  in  liquid  Helium.  The  solution  values  were  all 
obtained  from  RDFs  from  MC  whereas  the  X-state  parameters  have  been  verified  via  direct 
computation  of  moments  in  molecular  dynamics  (Eqs.  207-208)  as  well  as  using  RDFs  (Eqs.  209- 
210)  as  discussed  above.  Since  dynamical  simulations  for  the  He*  in  He  solution  are  difficult  due 
to  a  possible  reaction  between  He*  and  He,  such  a  comparison  (via  Eqs.  207-208)  is  not  available 
for  the  solution  at  the  present  time.  However,  the  pure  liquid  results  based  on  RDFs  give  credence 
to  the  solution  results  at  the  respective  pressures.  It  is  seen  that  the  parameters  are  very  sensitive  to 
the  electronic  states,  which  simply  reveals  that  the  primary  zone  must  ideally  include  more  than  the 
He*  atom.  The  table  also  contains  high  frequency  parameters  for  RDFs  from  each  given  electronic 
state  being  combined  with  potentials  for  different  electronic  states.  This  sort  of  comparison  reveals 
the  change  in  the  condensed  phase  dynamics  occurring  as  a  result  of  a  sudden  electronic  transition 
to  which  the  solvation  structure  has  not  adjusted. 

Figure  22  displays  a  comparison  of  pressure  dependent  VACFs  obtained  from  a  gaussian 
p(t)  model  with  MD-based  VACF  results  for  the  pure  liquid  helium.  The  gaussian  model  is  based 

2 

on  the  high  frequency  parameters  of  Table  5  and,  for  a  liquid  ( 00  e  =  (3(0),  since,  in  a  liquid,  one 
has  zero  adiabatic  frequency  for  a  tagged  particle),  is  given  by. 
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Figure  19.  Variation  of  fluid  structure  of  model  helium  with  pressure. 
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Figure  20.  Fluid  dynamics  (VACF)  of  model  helium  as  a  function  of  pressure: 

.  14.3  GPa 

.  12.7  GPa 

.  10.4  GPa 

O  1.4  GPa 
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Figure  21.  Variation  of  He(^S)  heatbath  modelling  (MTGLE)  parameters  with  pressure. 


Figure  22.  Comparison  of  pressure-dependent  He(^S)  VACF  obtained  from  MD  trajectories 
(continuous  line)  and  from  the  RDF-based  gaussian  friction  model  (dashed  line). 
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Figure  23.  Effect  of  pressure  on  He(^S)  VACF  from  RDF-based  gaussian  friction  model. 

_ 14  GPa 

10  GPa 

_  1.4  GPa 

Inset  shows  results  for  the  RDF-based  gaussian  friction  model  of  He(^S)  for 
comparison. 
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Table5a:  Some  Thermodynamic  Properties  of  liquid  4He  at  300K 


Vm  /  Cm3  mol _1 

(U>  /  kJ  mol*1 

(P)/GPa 

PV 

RT 

Bulk 

Solution 

Bulk 

Solution 

17.95 

0.36 

0.36 

0.29 

0.29 

2.2 

13.17 

0.74 

-2.0 

0.55 

0.55 

2.9 

0.85 

0.57 

9.02 

1.9 

-0.9 

1.4 

1.4 

5.1 

-2.2 

1.5 

4.76 

11.7 

10.4 

10.1 

10.1 

19.5 

11.9 

10.3 

4.35 

16.1 

14.2 

14.3 

13.8 

25.3 

15.6 

13.9 

(a)  In  the  case  of  solution  studies,  first  and  second  row  values  for  each  volume  correspond  to  a- 
and  c-states  of  He. 
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Table  6a:  GLE  parameters  for  He  and  He*  solutes  in  liquid  He 


P  (GPa) 

“eo  (THz) 

g>Ci(THz) 

X-state 

a-state 

c-state 

X-state 

a-state 

c-state 

0.55 

23,  (21) 

13 

23 

32,  (39) 

40 

27 

22 

15 

39 

19 

1.4 

32, (30) 

15 

32 

41,(48) 

35 

35 

32 

15  ■ 

38 

25 

10.1 

68, (65) 

___ 

48 

74,  (84) 

58 

53 

47 

— 

53 

54 

14.3 

77, (68) 

_ 

49 

83,  (86) 

62 

56 

48 

— 

55 

62 

(a)  The  MTGLE  parameter  values  for  He  are  given  as  direct  ensemble  averages  [see  Eqs.  (207), 

(208)  of  text]  with  the  values  in  parenthesis  being  from  the  corresponding  <p(r)  and  g(r)  [see  Eqs. 

(209) ,  (210)  of  text].  For  a-  and  c-states*of  He,  the  parameters  are  obtained  by  using  the 

corresponding  <p(r)  and  g(r)  data  and  are  reported  in  the  first  row.  The  second  row  values  for  a 
and  c  states  correspond  to  mixing  the  <p(r)  and  g(r)  functions  in  Eqs.  (209)-(210),  with  the  g(r) 

data  for  the  state  in  column  (a  and  c)  and  cp(r)  for  the  other  state  (c  and  a,  respectively)  to  which  the 
electronic  interaction  is  assumed  to  be  switched. 


P(t)  =(0g  exp 
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The  normalized  VACF,  denoted  X(t) ,  is  obtained  from  P(t)  by  solving  the  following  MTGLE 
response  function  equation: 

X(t)=-QJ20x(t)-Jp(t-t)%(x)dT  (2i2) 

where  dots  denote  time  derivatives.  The  results  shown  in  Fig.  22  are  based  on  the  Table  5 
parameters  obtained  via  the  RDF  route  (Eqs.  209-210).  It  is  seen  from  Fig.  22  that  the  gaussian 
model  improves  (by  agreeing  for  longer  times  within  the  short-time  region)  with  pressure,  though 
it  is  not  very  good  for  the  gaseous  phase  helium  corresponding  to  1.4GPa.  It  is  worse  at  0.5GPa 
(not  shown)  and  this  is  to  be  expected.  Although  the  gaussian  friction  model  based  on  the  direct 
moments  (cf.  Eqs.  207-208)  of  Table  5  led  to  overall  better  agreement  with  MD  VACFs,  but  we  do 
not  present  them  here.  The  gaussian  model  is  a  good  short  time  model  at  the  higher  densities  of  the 
liquid.  Anyway,  the  RDF  route  at  least  allows  us  to  obtain  the  VACFs  for  the  solution  of  He*  in 
He  from  the  RDF  data  for  He  bubble  reported  earlier.13  In  Figure  23,  we  show  the  VACFs  for 
the  solution,  based  on  the  gaussian  model  friction  for  c-state  helium  bubble  (obtained  using  the 
parameters  of  Table  5  and  Eqs.  21 1-212).  The  key  observations  are: 

(i)  the  excited  state  (c-bubble)  VACFs  qualitatively  change  with  increasing  pressure  in  a 
way  similar  to  the  ground  state  VACFs  (with  lOGPa  and  14GPa  having  essentially  a 
limiting  (saturated  effect)  behavior); 

(ii)  the  excited  state  VACFs  decay  on  longer  time  scales  than  the  ground  state  VACFs  of 
Fig.22  at  the  corresponding  pressures;  and 

(iii)  the  pressure  dependence  of  the  excited  state  VACFs  reflects  a  diminished  excited  state 
pressure  sensitivity  compared  to  the  ground  state  (see  Fig.  23  and  inset  for  comparison). 

(i)  is  a  natural  consequence  of  increased  density  of  the  liquid  leading  to  more  damping, 
whereas  (ii)  is  related  to  the  larger  size  of  the  excited  state  cage  (the  large  bubble  is  formed  due  to 
the  excited  electron's  repulsion  of  solvent  atoms),  (iii)  gives  a  significant  insight  into  what 
dynamical  effects  electronic  excitation  may  lead  to.  (iii)  is  due  to  the  farther  lying  atoms  of  the 
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liquid  playing  a  primary  role  in  the  solvation  effects  for  the  large  excited  state  atom's  (cage) 
bubble.  In  this  long-range  region,  RDFs  are  flat  and  less  sensitive  to  the  pressure,  and  hence  so 
are  the  moments  entering  the  gaussian  model.  This  trend  shows  that  long-time  properties  such  as 
diffusion  constants  are  probably  more  affected  by  pressure  dependence  of  condensed  phase 
influence  on  excited  states  than  the  present  short-time  dynamical  properties. 

4.  Conclusions 

The  present  pressure  dependent  simulations  provide  many  insights.  The  behavior  of  He  in 
the  condensed  phase  is  found  quite  different  from  the  case  of  Argon.  A  predominant  role  is  played 
by  entropy  effects  in  the  liquid  free  energy  of  binding  (reflected  in  the  the  internal  energy  of 
binding  being  an  excess  positive  quantity).  As  pressure  is  increased,  quantum  effects  become  less 
(eg.  at  lO-14GPa  cases  reported  here),  and  possibly,  they  may  become  relatively  important  only 
for  the  three  lower  pressures  (0.3- 1.4  GPa).  The  low  pressure  VACFs  indicate  gas-like  behavior 
with  a  long  hydrodynamic  positive  tail  whereas  the  high  pressure  ones  show  the  oscillatory 
structure  characterizing  a  liquid-like  condensed  phase  correlation.  The  RDFs  lead  to  similar 
conclusions. 

It  is  seen  that  the  MTGLE  parameters  are  quite  sensitive  to  (and  exponentially  dependent 
on)  the  pressure.  The  gaussian  friction  model  improves  with  pressure  and  provides  VACFs  for  the 
excited  atom's  bubble  from  structural  input  (RDFs)  alone.  The  sensitivity  of  the  dynamics  shown 
by  the  velocity  autocorrelation  functions,  for  He(3S)  and  He^S)  in  comparison,  reveals  the  change 
in  the  condensed  phase  dynamics  occurring  as  a  result  of  a  sudden  electronic  transition  to  which 
the  solvation  structure  has  not  adjusted.  Electronic  excitation  leads  to  the  sampling  of  lower 
frequency  components  of  the  solvating  medium  and  this  is  reflected  in  the  slower  decay  of  VACFs 
upon  electronic  excitation.  It  is  also  found  that  an  increase  of  pressure  may  introduce  less 
variations  in  the  translational  dynamics  of  the  excited  He(3S)  state  than  it  does  in  the  ground 
He('S)  state.  The  present  heatbath  models  are  being  employed  in  stochastic  dynamics  calculations 
of  electronically  inelastic  collision  processes21  within  the  high  pressure  4He  liquid. 
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Condensed  Phase  Dynamics 

Introduction 

The  condensed  phase  dynamical  problems  relevant  to  HEDM  research  involve  considerable 
richness  in  possible  complexity;  they  can  also  involve  quite  a  broad  range  of  phenomena  in  terms 
of  the  relevant  physics  and  chemistry.  Since  the  main  goal  of  the  three  year  program  was  to 
develop  new  and  required  methodologies  using  helium  metastables  as  prototypes,  our  specific 
calculations  have  been  determined  by  the  relevant  physics  of  this  system.  The  real  depth  of  the 
benefits  from  the  three  year  program  becomes  clear  only  by  imagining  the  intrinsic  generality  of  the 
computer  simulation  techniques  developed  here  and  the  fact  that  in  coordination  with  the 
simulations,  we  can  also  employ  previously  existing  technology  such  as  variational  transition  state 
theory1  to  treat  a  wide  variety  of  condensed  phase  reaction  and  aggregation  phenomena  that  have 
now  become  relevant  to  HEDM  preparation  and  storage  issues. 

The  helium  metastable  dynamics  in  a  high  pressure  helium  bath  can  be  adequately  described 
by  the  SCEM+GLE  technology  which  was  developed  by  this  effort  into  a  tool  for  chemistry  of 
species  embedded  in  the  bulk.  Previous  to  this  program,  the  existing  technology  was  mainly 
developed  for  gas-surface  collisions  and  assumed  asymptotic  noninteracting  species  undergoing  a 
collision  and  separating  again.  The  present  effort  began  by  completing  our  work  on  a  variable  time 
step  integration  of  the  stochastic  equations  of  motion  employed  in  the  SCEM+GLE  framework.2 

When  a  "collision"  occurs  within  the  bulk,  the  semiclassical  equations  do  not  lead  to  long 
time  asymptotic  constant  values  of  probabilities,  since  the  heatbath  is  constantly  interacting  with  the 
primary  "collision  partners"  and  perturbing  the  primary  zone  states  into  making  transitions.  We 
have  adapted  the  SCEM+GLE  methodology  for  liquid  bulk  by  (1)  implementing  a  thermal 
sampling  of  initial  conditions  for  the  "collision"  in  the  bulk  solution  and  (2)  developing  new  ways 
to  analyze  the  results  of  the  semiclassical  propagation  to  obtain  observables.  The  second  step 
presently  involves  examining  appropriate  time  correlation  functions.  In  particular,  probability 
fluctuation  autocorrelation  functions  yield  information  regarding  lifetimes  of  given  quantal  states 
from  their  decay  time.  Such  an  analysis  is  soundly  based  on  the  time  correlation  function 
formalism  of  nonequilibrium  statistical  mechanics.3 

The  next  subsection  describes  our  dynamical  calculations  using  the  SCEM+GLE 
technology  on  the  problem  of  metastable  helium  atoms  trapped  in  condensed  phase  bubbles.  It 
then  summarizes  the  overall  implications  of  our  results  in  terms  of  what  happens  to  the  stored 
energy  of  helium  metastable  atoms  in  gas  and  liquid  helium  as  a  function  of  pressure.  The 
subsection  also  finally  sketches  how  the  simulation  techniques  could  be  combined  with  other  in- 
house  techniques  at  Chemical  Dynamics  to  address  related  and  quite  different  possible  fates  of 
HEDM’s  in  condensed  phases;  this  type  of  framework  is  relevant  to  the  problem  of  Li  storage  in 
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hydrogen  or  other  conventional  fuels  such  as  hydrazine  which  are  being  presently  explored  at  the 
astronautics  laboratory  (AFSC). 

Dynamics  of  He*(3S)  in  Condensed  Phase  Bubbles 
1.  Solvated  Dynamics 

Since  the  metastable  helium  bubble,  whose  structure  was  discussed  in  the  previous  section, 
is  too  large  to  explicitly  include  all  the  several  dozen  solvent  He  atoms  of  even  the  first  solvation 
shell,  solvated  dynamics  was  modelled  by  one  solvent  on  either  side  of  a  collinear  triatomic  model, 
schematically  shown  in  Fig.  24,  where  the  central  atom  is  the  atomic  metastable,  He*(3S).  The 
dynamics  calculations  in  solution  were  based  on  employing  a  two  state  model  including  only  a  and 
X  states  along  with  their  spin-forbidden  radiative  coupling.  These  SCEM+GLE  simulations  were 
done  for  different  pressures  (in  the  GPa  range)  and  showed  no  significant  radiative  quenching  on 
time  scales  much  longer  than  that  found  for  reactive  clustering.  The  latter  timescale  was  estimated 
from  limited  MD  simulations  of  the  solution  system  whose  MC  study  is  reported  in  a  previous 
section  (see  p.  94).  The  former  was  manifested  by  a  nondecaying  probability  fluctuation 
autocorrelation  function  for  very  long  times.  The  net  conclusion  is  that  the  reactive  (clustering) 
channel  is  more  likely  than  radiative  quenching  in  the  bulk  and  hence  the  primary  loss  mechanism 
for  the  excited  atoms.  The  pressure  dependence  of  this  reactive  clustering  channel  was  revealed 
and  investigated  at  length  in  our  MC  simulations. 

STOCHASTIC  TREATMENT 
OF 

HELIUM  BUBBLE 


He  He*  He 


Figure  24.  Schematic  diagram  of  stochastic  model  with  an  elementary  primary  zone  made  of 
helium  metastable  atom  and  (effective)  "ghost"  atoms  that  represent  the  liquid 
helium  heatbath. 
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The  collinear  model  system  was  explored  numerically  at  this  stage  to  illustrate  the  basic  feasibility 
of  employing  the  present  methodology  and  to  display  the  appropriate  time  correlation  function 
behavior  expected  in  different  systems  with  different  condensed  phase  HEDM  lifetimes.  The 
collinear  triatomic  system  modelling  the  He*  bubble  in  a  high  pressure  He  liquid  matrix  was 
studied  using  the  semiclassical  eikonal  method  for  the  primary  zone  and  the  generalized  Langevin 
equation  (GLE)  to  simulate  heatbath  effects.  Probability  fluctuation  auto  correlation  functions 
(ACFs)  for  S  state  of  He  are  shown  in  Figs.  25  and  26  for  O.SGpa  and  1.4Gpa  a-bubble  which 
compare  well  to  the  displacement  autocorrelation  functions.  The  comparison  of  a  variety  of  mode 
ACFs  can  yield  information  about  the  roles  of  various  condensed  phase  modes  of  quenching  in 
destabilizing  the  HEDM.  The  generalized  Langevin  approach  correctly  incorporates  the  local  and 
macroscopic  mode  participations  in  primary  quenching  dynamics. 

The  dynamics  shown  is  for  radiative  quenching  of  He*(3S)  but  we  have  employed 
arbitrarily  increased  electronic  dipole  couplings  so  as  to  obtain  numerically  facile  rapid  quenching 
behavior.  As  mentioned  above,  no  radiative  quenching  is  observed  in  such  calculations  when  the 
actual  dipole  coupling  strength  of  Chabolowski  et.  al.4  is  employed  since  reaction  occurs  before 
quenching.  Special  dynamical  methods  are  being  developed  for  treating  quenching  in  very  long- 
lived  HEDM  species  (cf.  persisting  ACFs),  otherwise  requiring  long  trajectories.  The  technology 
employed  in  this  simple  model  may  also  be  applied  to  the  dynamics  of  vibrational  relaxation  in 
matrices;  such  a  semiclassical,  rather  than  classical,  description  is  especially  appropriate  if  V  — >  V 
transfer  is  an  important  relaxation  route.  The  adaptation  of  the  SCEM  scheme  can  be  accomplished 
by  expanding  in  a  suitable  joint  matrix-host  basis  set 

2.  Summary  of  Results  on  Helium  Metastables 

In  conclusion,  we  summarize  the  overall  picture  already  obtained  from  the  above  illustrative 
studies  of  the  helium  metastable  quenching  pathways  in  gas  and  condensed  phases. 

Radiative  versus  Nonradiative  Quenching  in  the  Gas  Phase 

Quenching  via  radiative  decay  in  vacuum  of  a  isolated  helium  atom  in  23S  state  is  controlled 
by  the  relativistically  induced  magnetic  multipole  couplings  within  the  atom.  Gas  phase  quenching 
at  higher  pressures  comes  from  He-He*  collisions.  Although  the  system  can  in  principle  quench 
nonradiatively,  it  was  found  that  the  related  probabilities  are  very  low,  ~10‘16,  compared  to  ~10-8 
for  the  radiative  routes.  The  direct  radiative  route  from  a  to  X  state  made  possible  by  spin- 
forbidden  radiative  coupling  seems  to  be  the  dominant  quenching  mechanism  in  all  our  studies  so 
far.  This  was  also  assumed  for  the  present  analysis  of  the  condensed  phase  problem. 
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He*(3S)  Probability 


P=0.5  Gpa  Displacement 


Figure  25.  Results  from  the  stochastic  dynamics  at  0.5  GPa  showing  the  fluctuation 
autocorrelation  function  (of  excited  He*);  this  function  will  display  decay  of 
excitation  due  to  electronically  inelastic  collisions  within  the  liquid  bulk,  (b)  the 
displacement  autocorrelation  function  of  the  metastable  He*  atom  in  liquid. 


He*(3S)  Probability 


P=1.4  Gpa  Displacement 


Time  (pS) 


Figure  26.  Same  as  Fig.  25  for  a  pressure  of  1.4  GPa. 

Radiative  Quenching  versus  Reaction  in  Solution 

The  behavior  of  an  excited  atomic  helium  species  in  condensed  phase  involves 
electronically  nonadiabatic  solute/solvent  interactions.  If  these  are  taken  to  be  reasonably  described 
by  the  gas-phase  two-body  potentials,  estimates  of  condensed  phase  behavior  may  be  obtained 
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from  the  computer  simulations  described  in  the  previous  subsections.  The  extremely  small 
nonradiadve  and  many-state  collisional  quenching  pathways  of  Fig.  4  are  assumed  to  be  ignorable 
upon  transfering  to  condensed  phase  in  our  present  analysis.  This  means  two-state  models,  such 
as  the  a  and  X  combination  used  in  Fig.  5,  are  adequate  in  the  condensed  phase  and  hence  the  a 
and  c  states  are  the  only  key  potentials.  Interactions  of  a  He*(2^S)  in  a  high  pressure  helium  liquid 
lead  to  quantitatively  different  stabilities  for  the  metastable  atomic  species  depending  on  the  a-state 
and  c-state  characters  of  solute-solvent  potentials.  The  c-state  bubble  protects  the  atomic  species 
up  to  very  high  pressures  whereas  the  a-state  more  easily  allows  reaction  with  a  solvent  He  atom 
leading  to  a  pronounced  degree  of  cluster  formation  that  varies  with  pressure. 

One  remarkable  feature  (seen  in  the  Monte  Carlo  simulation  of  helium  bubble  states)  is  the 
maintenance  of  a  bubble  structure  despite  the  cluster  formation  (at  least  when  the  simulation 
employs  reasonable  pairwise  additive  interactions).  The  c-  and  a-states  are  radiatively  coupled  in 
reality  and  both  states  can  play  a  role  in  determining  the  metastable  atom's  equilibrium  bubble 
shape.  The  cluster  formation  observed  here  is  a  neutral  excitation  trapping  process  in  liquid  that  is 
analogous  to  the  charged  exciton  trapping  phenomenon  observed5  in  rare  gas  halide  crystals.  These 
results  indicate  that  stabilization  of  a  metastable  by  forming  clusters  of  it  with  a  suitable  co-species 
may  still  result  in  trapping  into  bubbles  in  the  condensed  phase.  Since  such  bubbles  can  protect  the 
species'  energy  from  dissipation  to  bulk,  it  is  useful  to  identify  conditions  for  their  formation  in 
different  matrices  for  different  species. 

It  is  verified  that  reaction  is  the  dominant  loss  mechanism  in  condensed  phases  when 
compared  to  collision-induced  radiative  decay.  It  was  found  that  the  reaction  can  go  on  to  result  in 
formation  of  higher  clusters,  but  sensitively  depending  on  pressure  and  electronic  state  of 
interactions.  The  solution  structure  has  been  examined  via  computer  simulations.  The 
parameterization  for  a  GLE  heatbath  description  obtained  that  were  then  used  in  exploratory 
condensed  phase  dynamics  of  the  HEDM.  It  has  been  found  that  there  is  no  radiative  quenching  on 
the  time  scales  within  which  significant  reaction  occurs. 

Cluster  Stability  and  Concentration  Effects  in  Solution 

Clearly,  the  above  results  show  that  the  average  degree  of  aggregation  achieved  by  the 
clustering  events  depends  on  the  thermodynamic  state.  Another  factor  in  controlling  their  stability 
is  the  concentration  of  such  clusters  each  of  which  could  potentially  release  the  stored  energy  in 
collisions.  These  studies  require  many-body  potentials  and  we  do  not  have  any  information 
pertaining  to  these  yet.  But  it  is  believed  from  concentration  dependence  studies  that  dimer-dimer 
collisions  in  liquid  reduce  dimer  lifetimes  to  as  small  as  30ms6  at  a  density  of  1012cm*3. 
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Extended  Studies 


Solvent  Shift  of  2^S->2^P  Absorption  Line 

Another  key  application  of  the  present  technology  is  to  obtain  spectral  absorption  and 
emission  profiles  for  HEDM's  in  the  condensed  phase  that  include  detailed  solvent  static  and 
dynamical  shifts  and  any  shape  changes.  Such  calculations  can  be  valuable  for  interpreting  the 
mechanisms  controlling/modifying  condensed  phase  experimental  spectral  data  and  in  comparing  to 
gas  phase  spectra.  We  have  performed  illustrative  calculations  using  the  recent  helium  potentials  of 
Yarkony  to  reveal  this  capability. 

One  of  the  early  key  experimental  observations  on  helium  bubbles  was  the  blue  shifted 
nature  of  the  above  absorption  line  in  contrast  to  a  red  shifted  nature  of  the  corresponding  emission 
line.  This  was  fruitfully  studied  in  a  phenomenological  bubble  model7  using  adiabatic  line  shape 
theory  in  the  static  limit.  Solvent  static  and  dynamical  effects  can  both  in  general  play  a  role  in 
shifting  spectra.  Although  we  have  characterized  the  high  pressure  helium  bath  and  deduced 
parameters  for  stochastic  dynamics  simulations  of  metastable  species  trapped  in  it,  in  the  case  of 
the  present  atomic  transition,  solvent  dynamics  can  be  shown  to  be  less  crucial  in  determining  the 
absorption  line  at  high  pressures  (a  typical  bath  atom  moves  only  ~0.2aQ  within  a  transition  dipole 
relaxation  time  compared  to  ~1.5a0  for  the  low  pressure  cryogenic  regime  where  previous 
experiments  were).  We  have  verified  from  our  studies  the  previous  experimental  absorption  line 
blue  shifts  employing  the  quantal  radial  distribution  function  from  Hansen  and  Pollock8  and 
Yarkony's1  transition  dipole  function.  We  have  also  performed  a  similar  calculation  employing 
our  radial  distribution  function  for  the  a-bubble  and  the  a<->b  transition  dipole  function  to  find  a 
large  red  shift  for  the  high  pressure  line  position  by  ~.29eV.  This  value  should  be  modified  by 
considering  contributions  due  to  transitions  dictated  by  the  (presently  unknown  radiative)  coupling 
of  the  c-state  to  the  u-manifold  (outside  the  primary  space  employed  by  Yarkony)  of  the  25P 
asymptote.  We  have  found  that  the  large  red  shift  is  predicted  even  when  the  contribution  from  the 
clustering  region  of  a-bubble  is  subtracted.  At  slightly  lower  pressures,  the  complication  of 
reaction  may  be  avoided  and  the  degree  of  a,c  mixing  characterized  via  spectrosopic  studies. 

Studies  on  the  H2  matrix 

We  have  conducted  preliminary  studies  on  condensed  phase  liquid  hydrogen  (high 
pressure,  in  the  Gpa  range,  and  hence  classically  described)  with  a  two  center  Lennard- Jones 
(2LJC)  anisotropic  effective  pair  potential.  The  goal  of  using  an  anisotropic  potential  even  though 
liquid  hydrogen  is  known  to  be  quite  isotropic  in  behavior  was  to  develop  the  codes  and  analysis 
tools  for  such  a  general  diatomic  condensed  phase  matrix. 
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The  first  step  is  to  devise  an  effective  anisotropic  pair  potential  for  condensed  phase 
simulations  of  molecular  hydrogen.  MD  calculations  have  been  earned  out  for  liquid  H2  system  of 

108  molecules  interacting  through  two  Lennard- Jones  centers  (2UC)  coincident  with  the  positions 
of  the  atomic  masses.  Figures  27-29  display  the  center  of  mass  (c.m.)  pair  correlation  functions 
(PCFs)  for  special  configurations  and  the  time  correlation  functions  (TCFs)  for  2LJC-model  H2  at 

265K  and  4  GPa. 


Liquid  Hydrogen  (2LJC) 


Figure  27.  Simulation  results  for  radial 

distribution  function  of  parallel  (P) 
and  L-shaped  hydrogen  dimers 
in  liquid. 


Liquid  Hydrogen  (2LJC) 


Figure  28.  Simulation  results  for  radial 
distribution  function  of  T- 
and  X-  shaped  hydrogen  dimers 
in  liquid. 


Figures  27-28  display  c.m.  PCFs  for  configurations  which  lie  within  ±  10^  of  specific  relative 
orientations:  ’X'='crossed',  'P'='parallel  adjacent',  T=’T-geometry’,  and  ,L’’="linear,  or  parallel 
end  to  end’.  The  position  of  the  maxima  of  these  PCFs  compare  well  with  the  position  of  the 
potential  energy  minima  for  two  isolated  molecules  in  the  corresponding  configuration.  This 
indicates  that  the  presence  of  other  molecules  in  the  dense  liquid  have  little  effect  on  the  minimum 
of  the  potential  field  which  acts  on  adjacent  molecules.  By  contrast,  for  the  models  of  other  related 
linear  molecules,  (K.  Singer,  A.  Taylor,  and  J.  V.  L.  Singer  Mol.  Phys.  33,  1757  (1977))  the 
same  arrangements  have  been  found  to  be  almost  equally  stable.  Figure  29  displays  the  c.m.  and 
angular  velocity  auto-correlation  functions  of  the  2UC-model  H2.  In  the  liquid  state,  the  VACF 

exhibits  a  negative  minimum,  which  is  interpreted  as  caused  by  bac  . -scattering  by  nearest- 
neighbors,  and  a  long  negative  tail,  ascribed  to  a  cooperative  motion  of  the  surrounding  particles. 
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The  double  minimum  in  the  c.m.  VACF  can  be  interpreted  (J.  Barojas,  D.  Levesque,  and  B. 
Quentrec  Phys.  Rev.  A2,  1092  (1973))  as  arising  from  the  different  mobility  of  a  linear  molecule 
parallel  and  at  right  angles  to  its  axis.  This  is  clear  from  the  c.m.  VACFs  of  velocity  components 
parallel  and  perpendicular  to  the  molecular  axis.  The  minimum  of  VACF  is  deeper  and  occurs  at 
earlier  times  than  that  of  VACF||.  The  presence  of  double  minimum  in  VACF  was  interpreted  to 

be  the  result  of  successive  negative  impulses  transmitted  to  the  c.m.  first  by  one  and  then  by  the 
other  end  of  the  molecule,  as  it  encounters  obstacles.  The  angular  VACF  would  also  exhibit 
negative  minimum  for  more  anisotropic  systems.9 


Liquid  Hydrogen  (2LJC) 


Time  (pS) 

Figure  29.  Simulation  results  for  autocorrelation  functions  of  various  components  of 
the  molecular  velocity  in  liquid  hydrogen;  correlation  functions  of  center  of 
mass,  parallel,  perpendicular  and  angular  velocities  are  shown. 

Detailed  Studies  of  Molecular  Structure: 

The  pair  correlation  function  (pcf)  for  linear  molecules  is  generally  defined  as 
g(R,  9  j,  0  ^  9  j2) ,  where  ^  =  12I  ’  ^  12  vector  joining  the  two  centers  of  mass; 

A  A  A  A 

cos  0  J  =  1 J  •  R 12,  cos  0  2  =  1  2  •  R 12,  and  cos  <p12  =  (I  x  R  2)  •  (1  2  x  R  12>;  I  r  I  2  are 

the  bond  vectors  of  molecules  1  and  2,  and  A  denotes  unit  vector.  It  is  difficult  to  study  the 
properties  of  such  a  4-dimensional  correlation  function.  The  expansion  in  terms  of  spherical 
harmonics  can  be  used  to  analyze  the  pcf  of  linear  molecules10; 
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(213) 


g(8,  8  r  9  2,  <p12)  =  4*  ZSIsU  m  «  \  m<«r  «>PY1.  -  m<«?  f  2> 


where  Y  ^  m>  Y  _  m  are  normalized  associated  Legendre  functions.  Although  not  directly 
measurable  until  very  recently,11  the  coefficients  gj  j  m  can  be  obtained  from  the  translational 
and  rotational  trajectories  of  the  molecules  in  the  computer  "experiment". 


The  site  (atom)-site  pcf  g  s(r)  defined  by 

,2 


^T7  g  s(  r )  =  lim  y  SSZSn  (r  <r.  .  <  r  +  Ar )  / 47tr 2  Ar 

2V  Ar  ->0  1  i  s  j  s  1SJS 


(214) 


where  N  is  the  number  of  molecules,  V  is  the  volume,  s  and  s'  correspond  to  the  nonbonding  sites 
of  the  pair  of  molecules  i  and  j,  and  n(r)  is  the  histogram  count.  gs(r)  is  related  by  Fourier 
inversion  to  the  static  structure  factor  which  can  be  determined  by  X-ray  or  neutron  diffraction. 
Useful  information  has  also  been  obtained  from  the  incidence  of  certain  special  configurations  as 
functions  of  molecular  separation. 

The  site-site  pcf,  the  first  six  coefficients  in  the  above  expansion  and  the  four  special 
configurations  have  been  calculated  for  molecular  H2  at  270K  (~4  GPa).  The  principal 

characteristics  of  these  structural  features  are  reported  in  Table  7. 

Curves  for  the  first  six  radial  distributions  which  arise  as  non-vanishing  coefficients  in  the 

expansion  (1)  are  shown  in  Figs.  30  and  31.  The  general  appearance  of  the  g  j  j  m is  surprisingly 
similar  to  the  long  (more  anisotropic)  linear  molecules.12. 

gooo(R)  1S  the  c.m.  pcf.  The  main  peak  occurs  at  5.83  bohr  which  coincides  with  the 
peak  of  the  corresponding  T  pcf.  The  second  peak  is  at  somewhat  less  than  twice  the  distance  of 
the  first  peak,  Rmax.  It  could  obviously  arise  from  two  vectors  of  length  Rmax  making  a  small 
angle  with  each  other.  There  is  no  bulge  on  the  left-hand  side  of  the  first  peak.  This  is  mostly  the 
case  with  less  anisotropic  molecules  such  as  N2,  F2.  Generally  its  range  of  correlation  has  been 

found  to  be  greater  than  that  of  gs(r).  As  we  employed  108-molecules  in  our  preliminary 
simulation  studies,  the  correlation  length  is  restricted  and  does  not  extend  beyond  the  second 
maximum.  The  necessity  to  carry  out  large  system  is  thus  apparent 
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Table  7:  Characteristics  of  pair  correlation  functions 
Positions  (bohrs)  of  maxima  of  pcf  s  for  special  configurations 

X  P  T  L 

5.29  5.41  5.83  6.43 


Maxima  at 

Minima  at 

(bohr) 

(bohr) 

gs(r) 

5.69,  10.57 

8.09 

g000 

5.83,  10.6 

8.05 

g  200 

6.28 

5.41 

g220 

5.23,  6.67 

5.71 

g2  21 

6.19 

5.29 

®200  *  cos^  ®  i  “  *)  has  no  negative  lobe  at  small  R;  this  is  because  if  the 

R 

centers  of  mass  are  close  to  one  another  the  vector  connecting  them  is  likely  to  be  approximately 
normal  to  the  axes  of  both  molecules  -as  in  the  configurations  X  and  P.  The  maxima  of  these 
pcf  s  is  in  fact  close  to  the  minimum  of  g  ^ qq-  The  zeros  of  g  £qq  result  from  the  cancellation  of 


positive  and  negative  contributions  to  ^3  cos^  9  ^  -  l)  arising  from  0  ^  >  cos  ^(a/3  /3 ); 

they  occur  at  separations  which  are  close  to  the  main  peaks  of  Sqqo-  For  the  T  configuration 

1 3cos^0j  -  1^  equals  either  2  or  -1,  and  a  small  mean  value  of  the  angular  factor  for 

arrangements  close  to  T  is  plausible.  This  is  compatible  with  a  large  contribution  of  T 
configurations  to  the  main  peak  of  Sqoo- 


145 


Figure  30.  Simulation  results  for  gooo  §200*  ^  §220  coefficients  of  Eq.  (213). 
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g«'m  (R) 


gN.m  coefficients 


R/bohr 
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Figure  31.  Simulation  results  for  §221*  8222 >  8400  coefficients  of  Eq.  (213). 
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For  g  22Q  *  f  (R)^3  cos^  0  ^  -  1^3  cos^  0  ^  -  l)  ,  the  initial  positive  lobe  with  a 

R 


maximum  at  5.23  bohr  indicates  0  ^  =  0  2  =  ft  /  2 ,  i.e.  +  and  II  like  configurations.  The  lowest 

possible  value  of  the  angular  factor  occurs  at  0  ^  =  tc  /  2,0  ^  =  0  or  vice  versa,  i.e.  for  T 

configuration.  Table  7  shows  that  the  minima  of  g  ^2q  almost  coincides  with  the  maxima  of 

§  000-  This  again  points  to  a  significant  incidence  of  T  configurations  at  the  main  peak  of  §  QOO . 
The  subsequent  positive  lobe  is  probably  due  to  L-like  configurations. 

The  geometrical  interpretation  of  the  higher  expansion  coefficients  becomes  rapidly  more 
complicated.  The  theoretical  significance  of  the  expansion  lies  in  its  possible  use  in  perturbation 
expansion  for  anisotropic  potentials. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


From  our  gas  phase  method  development  efforts,  based  on  the  semiclassical  computational 
approach  (within  the  eikonal  approximation)  we  conclude  that  the  approach  is  useful  for 
applications  to  collisional  problems  that  are  difficult  to  treat  quantum  mechanically  due  to  the  large 
number  of  nuclear  rovibrational  states  (coupled  channels).  The  present  semiclassical  methodology 
reveals  much  formal  flexibility,  and  techniques  are  now  available  for  applying  it  to  different  levels 
of  treatment  of  both  classical  and  quantum  dynamics  of  polyatomic  systems.  For  each  application, 
the  principles  validated  by  our  research  should  be  carefully  exploited  in  choosing  a  computational 
strategy;  the  strategy  should  depend  on  the  nature  of  the  potential  energy  surface  on  which  the 
collision  dynamics  evolves,  the  accuracy  of  the  available  potentials  and  couplings  and  the  accuracy 
with  which  state-to-state  results  are  needed.  Validation  studies  described  in  this  report  are  not 
exhaustive,  but  do  illustrate  of  the  potential  of  the  methods.  New  directions  for  novel  and  efficient 
semiclassical  computations  have  been  formally  developed  by  employing  semiclassical 
wavefunctions  for  multidimensional  electronically  inelastic  problems,  although  the  validations  of 
this  new  approach  (studies  of  O+HF  and  Na+H2  collisions)  have  been  done  only  for 
multidimensional  single  potential  surface  problems.  A  "primitive"  version  of  the  theory  (not  using 
semiclassical  wavefunctions)  was  employed  for  most  of  this  research,  since  the  formal 
improvements  were  discovered  towards  the  end  of  the  research.  The  primitive  method  has  been 
illustrated  using  the  (many-surface)  electronically  inelastic  problem  of  metastable  helium  collisions 
in  both  gas  and  condensed  phases. 

The  effort  on  condensed  phase  modelling,  and  simulations  on  helium  metastable  structure 
and  dynamics  in  liquids  illustrate!  the  power  of  the  computational  approach  by  revealing 
microscopic  aspects  of  the  phenomenon  of  bubble  formation  around  an  electronically  excited  atom 
in  a  liquid.  The  sensitivity  of  the  simulation  method  to  the  potential  surfaces  for  different  electronic 
states  shows  that  semiquantitative  levels  of  treatment  are  possible.  Even  using  approximate 
potential  energy  surfaces,  such  calculations  can  be  used  to  reveal  mechanistic  details  that  are  not 
easily  accessed  in  experiments,  and  assist  in  the  interpretation  of  condensed  phase  chemistry 
experiments.  Our  experience  in  building  a  stochastic  model  for  the  simple  helium  system  can  lead 
the  way  for  treating  other  inert  gas  hosts;  however,  heavier  inert  gas  systems  do  exhibit  different 
phenomena  and  require  more  complex  potential  surface  information.  Systems  such  as  the  hydrogen 
matrix,  also  recently  studied  by  us,  can  display  more  complex  behavior  that  may  be  difficult  to 
obtain  accurate  potential  energy  information  on;  this  can  happen  if  the  excited  state  chemistry 
involves  many  hydrogen  molecule  neighbors,  and  atom  transfer.  The  key  requirement  for 
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dynamical  studies  is  potential  energy  surface  information,  and  systems  for  which  this  is  available 
are  amenable  to  theoretical  studies. 

The  ground  state  dynamics  of  different  species  in  the  hydrogen  matrix  can  be  treated  by 
combining  the  present  methods  with  more  standard  methods  such  as  variational  transition  state 
theory.  Besides  computing  rates  for  elementary  chemical  processes  and  lifetimes  of  species,  the 
technology  developed  here  is  also  useful  for  obtaining  theoretical  line  shapes  for  solid  matrices,  a 
quantity  directly  probed  in  experiments. 
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Byron  H.  Lengsfield  ID  and  George  Adams'  research  group 

There  were  close  collaborations  between  Chemical  Dynamics  and  the  Ballistics  Research 
Laboratory  (BRL)  during  this  period.  Several  visits  and  discussions  occurred.  BRL  calcula!:ons, 
in  collaboration  with  the  Johns  Hopkins  effort  of  David  Yarkony,  provided  information  on  the 
potential  energy  surfaces  for  helium  and  other  helium  and  hydrogen  metastable  species.  The  He-He 
system  was  discussed  in  depth. 

David  R.  Yarkony 

There  were  fruitful  interactions  on  the  helium  problem  and  the  requirements  expected  of 
quantum  chemistry  for  solving  electronically  inelastic  radiative  and  nonradiative  dynamical 
problems. 

Daniel  Konowalow 

The  main  interactions  were  on  the  sensitivity  of  condensed  phase  effects  on  long  range 
helium  potentials  and  the  possibility  of  obtaining  ab  initio  potentials  for  use  in  L1/H2  matrix 
simulation  studies. 

Marcy  Rosenkrantz 

The  main  interactions  were  on  methods  for  different  species-matrix  interactions  to  be 
employed  for  dynamics  studies. 

Steve  Rodgers 

A  number  of  significant  interactions  occurred  regarding  the  important  technologies  for 
HEDM  research,  especially  for  condensed  phase  dynamical  problems. 
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Ara  Apkarian 

Stimulating  discussions  took  place  including  a  visit  by  Dr.  Apkarian  to  Chemical  Dynamics 
to  discuss  possible  theoretical  support  of  his  experimental  program. 

Nathan  Presser 

The  main  interaction  was  on  the  condensed-phase  line  shape  problem  of  sodium  in  rare  gas 
solids,  a  problem  to  which  the  technology  developed  on  this  project  can  be  applied. 

Gilbert  Collins 

Recent  interactions  have  occurred  on  the  problem  of  hydrogen  atom  diffusion  and  reaction 
in  hydrogen  matrices  and  how  to  interpret  the  experiments  on  thermal  spikes  and  isotope  effects. 

The  following  Chemical  Dynamics  staff  members  have  been  partially  supported  during  this 
three  year  contract: 

Dr.  Pazhayannur  K.  Swaminathan 

Dr.  Michael  J.  Redmon 

Dr.  Cheruvu  S.  Murthy 

Dr.  Bruce  C.  Garrett 
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